Staticity, Self-Similarity and Critical Phenomena 
in a Self- Gravitating Nonlinear a Model 



Dissertation 
zur Erlangung des akademischen Grades 
Doktorin der Naturwissenschaften 
an der Formal- und Naturwissenschaftlichen Fakultat 
der Universitat Wien 



eingereicht von 
Christiane Lechner 



im Oktober 2001 



Kunst ist schon, 
macht aber viel Arbeit. 

Karl Valentin 



ZUSAMMENFASSUNG 



Ziel dieser Arbeit ist eine dataillierte Untersuchung spezieller Aspekte der Dy- 
namik einer Klasse nichtlinearer Materiefelder im Rahmen der allgemeinen Rel- 
ativitatstheorie. Von besonderem Interesse hier sind statisdie Losungen (unter 
Beriicksichtigung einer positiven kosmologischen Konstante), selbstahnliche Lo- 
sungen sowie die Bildung schwarzer Locher. 

Selbstgravitierende Materiefelder zeigen fiir Anfangsdaten, die an der Schwelle 
zum Kollaps liegen, ein Verhalten (sogenannte "kritische Phanomene"), das durch 
"Scaling", Selbstahnlichkeit (bzw. Statizitat) und Universalitat charakterisiert 
ist . Das hier - in spharischer Symmetrie - untersuchte SU(2) a Modell, ist in 
diesem Zusammenhang von besonderem Interesse, da die dimensionslose Kop- 
plungskonstante in nichttrivialer Weise in die Theorie eingeht. Ziel dieser Arbeit 
ist, kritische Phanomene, insbcsondere den Skalcncxponenten und die kritische 
Losung, in Abhangigkeit der Kopplung zu untcrsuchcn. 

Die Untersuchung erfolgt in zwei Schritten. Zunachst werden selbstahnliche 
Losungen (unter Zuhilfenahme der Symmetrie) numerisch konstruiert und deren 
Stabihtat untersucht (Kapitel 4). Wir reproduzieren die Ergebnisse von Bizon 
et al., die fiir kleine Kopplungen eine einparametrige Familie von kontinuierlich 
selbstahnlichen (CSS) Losungen konstruiert haben. Wir zeigen, dafi die erste 
Anregung dieser Familie eine instabile Mode hat. Welters konstruieren wir eine 
diskret selbstahnliche (DSS) Losung fiir grofie Kopplungen, die ebenfalls eine in- 
stabile Mode hat. Wir stellen die Hypothese auf, daB die DSS Losung bei einem 
bestimmtcn Wert der Kopplung aus der ersten CSS Anregung in einer homoklinen 
Loopbifurkation entsteht. 

Im zweiten Schritt werden einparametrige Familien von Anfangsdaten numerisch 
in der Zeit zu entwickelt (Kapitel 5). Nahekritische Anfangsdaten erhalt man 
durch Bisektion. Die kritische Losung kann an dem (zeitlich) intermediaren Ver- 
halten der nahekritischen Evolutionen abgelesen werden. Der Skalenexponent 
ist durch die Masse der schwarzen Locher in Abhangigkeit des Parameters der 
Anfangsdaten bestimmt. Die so erhaltenen Ergebnisse stimmen sehr gut mit 
den Eigcnschaften der oben konstruicrtcn selbstahnlichen Losungen iiberein. Fiir 
kleine Kopplungen ist die kritische Losung die erste CSS Anregung. Fiir groBe 
Kopplungen ist es die DSS Losung. Fiir mittlere Kopplungen finden wir einen 
Ubergang von CSS zu DSS in der kritischcn Losung. Dieser Ubergang ist mit der 
Hypothese der homoklinen Loopbifurkation konsistent. 

Zusatzlich wird in dieser Arbeit iiber statische Losungen des Modells (mit posi- 
tiver kosmologischer Konstante) berichtet (Kapitel 3). 
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Abstract 



The aim of this work is to study certain aspects of the dynamics of a class of 
self-gravitating non-linear matter fields. In particular we concentrate on soli- 
ton solutions (in the presence of a positive cosmological constant), self-similar 
solutions and the formation of black holes. 

The dynamics of initial data at the threshold of black hole formation are char- 
acterized by phenomena (so-called critical phenomena) including scaling, self- 
similarity (resp. staticity) and universality. In this work we concentrate on SU(2) 
a models coupled to gravity in spherical symmetry. These models are interest- 
ing due to a dimensionless parameter - the coupling - which enters the theory 
non-trivially. The aim is to investigate how critical phenomena ~ in particular 
the critical solution and the scaling exponent - depend on the coupling. 

We use two essentially different methods to study the threshold behavior: Making 
use of the symmetry both discrete (DSS) and continuous (CSS) self-similar solu- 
tions are constructed (numerically) by solving boundary value problems (Chapter 
4). The stabihty of these solutions is studied. For small couphngs, reproducing 
results of Bizon et al., we find a discrete one-parameter family of CSS solutions. 
Of particular interest is the first CSS excitation. For large couplings we construct 
a DSS solution. Both solutions have one unstable mode. We conjecture, that the 
DSS solution bifurcates from the CSS solution in a homoclinic loop bifurcation 
at some value of the coupling constant. 

The second method consists of evolving one parameter families of initial data 
numerically (Chapter 5). By a bisection search the initial data are fine-tuned 
such that they are close to the threshold. The critical solution then is determined 
by the intermediate asymptotics of near-critical data. The scaling exponent is 
determined from the black hole mass as a function of the parameter in the initial 
data. Our results are in very good agreement with the results on the self-similar 
solutions we obtained above. For small couplings the critical solution is CSS, for 
large couplings it is DSS and for intermediate couplings we find a new transition 
from CSS to DSS in the critical solution, which shows "episodic CSS". This 
transition is consistent with the hypothesis of the homoclinic loop bifurcation. 

In addition this work also contains results on static solutions of the model in the 
presence of a positive cosmological constant A (Chapter 3). 
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Chapter 1 
Introduction 



The theory of general relativity describes gravity in terms of the curvature of 
the four dimensional Lorentzian manifold representing spacetime. The Einstein 
equations, G^^, = i^T^y, relating matter to the curvature of spacetime involve the 
geometric objects T^jy, the stress energy tensor of matter, and G^y, the Einstein 
tensor. Behind this simple geometric formulation there hides a coupled system 
of ten quasi-linear partial differential equations (PDEs), with the components of 
the metric as dependent variables. In addition, the matter fields themselves are 
subject to some field equations, which in turn contain the metric and its deriva- 
tives. These field equations therefore have to be solved simultaneously with the 
Einstein equations. Of special physical interest is to formulate and solve the ini- 
tial value problem. In the "3-1-1" formulation for example, initial data consist 
of two symmetric three-tensor fields given on an initial spacelike hypersurface, 
which are subject to four constraint equations. The remaining six Einstein equa- 
tions are used to evolve these initial data in time. One question of interest is for 
example: given smooth initial data, what is the long time behavior of the solu- 
tion? There are analytic results showing global existence for sufficiently "weak" 
initial data^. On the other hand, singularity theorems^ predict that sufficiently 
"strong" initial data develop a singularity in finite time. If the cosmic censorship 
hypothesis holds, these singularities should be shielded by a horizon such that 
they are invisible to distant observers. However, according to the complexity of 
the equations it is clear that i) only a very small number of exact solutions is 
known and ii) it is very difficult to get an analytic handle on the equations. It is 
therefore both necessary and extremely fruitful to combine analytic approaches 
with a numerical treatment of the equations. 

An important step to understand the dynamics is to study the possible "end 
states" for the given matter model. It is reasonable to assume that regular initial 
data will asymptote to a stationary solution at late times, which could be e.g. a 

^In particular For a review see |65| . 
^See inn and [SI. 
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stationary stable black hole, a stable soliton solution (corresponding to a star) or 
dispersion leading to Minkowski spacetime. Stationarity reduces the equations to 
an elliptic system. If additional symmetries are imposed, the problem simplifies 
further. In particular a static soliton or black hole solution, which is spherically 
symmetric, is obtained by solving a coupled system of (nonlinear) ordinary differ- 
ential equations (ODEs). Initiated by the work of Bartnik and McKinnon p], who 
numerically construced static soliton solutions to the Einstein Yang-Mills (EYM) 
system in spherical symmetry, several self-gravitating matter models have been 
studied with regard to static soliton or black hole solutions^. Unfortunately most 
of the (non-trivial) solutions found are unstable, so they are not relevant for the 
late time behavior of general initial data. However, a few years later, it turned 
out that static solutions with a single unstable mode play an important role as 
intermediate attractors in type I critical collapse, as explained below. 

Another fascinating field of research was started by the work of Choptuik P^ ll9j. 
who investigated the threshold of black hole formation for the self-gravitating 
massless Klein-Gordon system in spherical symmetry. Due to the analytic work 
by Christodoulou^ it was known that "small" (in a well defined sense) initial data 
disperse, such that in the long time evolution such data approach Minkowski 
spacetime, whereas "strong" initial data form a black hole. Choptuik numeri- 
cally evolved one parameter families of initial data, that interpolated between 
black hole formation and dispersion. He fine-tuned the parameter such that a 
"tiny" change in this parameter changed the end state from black hole formation 
to dispersion or vice versa . Speaking in the language of dynamical systems, the 
space of initial data (of the model under investigation) is divided into basins of 
attraction, the attractors in this case being black holes and Minkowski space- 
time. From this point of view, Choptuik studied the boundary between two such 
basins of attraction. The phenomena he found, called critical phenomenal , can 
be summarized by the keywords scaling, self- similarity and universality. Scaling 
relates the black hole mass to the parameter p in the initial data via the simple 
power law niBH oc (p — p*y, where p* is the critical parameter of the family. In 
particular this means that the black hole masses can be made arbitrarily small 
by fine-tuning the initial data. Furthermore, all near critical evolutions approach 
a self-similar solution at intermediate times. Studying several families of initial 
data in this way Choptuik found that these phenomena, especially the scaling 
exponent 7 and the self-similar solution are independent of the family. So critical 
phenomena are universal within a given model. The numerical resolution of the 
phenomenon of self-similarity required sophisticated methods. Encouraged by 

■^Most of the work was done in spherical symmetry. The Einstein Yang-Mills system was also 
investigated in axis symmetry. In addition there are investigations concerning slowly rotating, 
that is stationary solutions to this model, obtained as linear rotational perturbations of the 
Bartnik McKinnon solutions and the coloured black holes. For details and an extensive list of 
references see the review article by Volkov and Gal'tsov |71] 

*For references see e.g. |57] . 

^For details and review articles see Chapter |S| 
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Choptuik's results several other models were investigated with respect to critical 
phenomena. Most of the work concentrated on spherical symmetry with the early 
exception of Abrahams and Evans P] , who considered axially symmetric gravi- 
ational waves. All models exhibit the phenomena described above. The scaling 
exponent 7 and the self-similar solution turned out to depend on the model. In 
particular there are models were the intermediate attractor is continuously self- 
similar (CSS) and other models, where it is discretely self-similar (DSS). Up to 
now it is not clear what causes the symmetry to be discrete or continuous. The 
results by Bizon et al. ^U], who studied threshold phenomena of the SU(2) a 
model in flat space, indicate that critical phenomena are not tied to the Einstein 
equations, but are rather a general feature of hyperbolic PDEs. On the other 
hand, up to now it is not clear whether there are models in flat space which allow 
for a discretely self-similar solution, or whether DSS is a special feature of the 
Einstein equations. 

It is worth noting that Choptuik's numerical work contributed much to the gen- 
eral understanding of the Einstein equations. In particular it stimulated the 
(semi-analytic) study of self-similar solutions and their stability properties. It 
turned out that the self-similar solutions at the threshold of black hole formation 
have one unstable mode. Furthermore a scaling argument relates the eigenvalue 
A of the unstable mode to the scaling exponent 7 via 7 = 1/A. Taking these 
results together, one has a good (although not rigorous) understanding of how 
critical phenomena emerge. 

Apart from the phenomena found by Choptuik, some models in addition give rise 
to another type of critical behavior. There the intermediate attractor is a static 
(or oscillating) solution with one unstable mode. The black hole masses formed by 
slightly super-critical data are not arbitrarily small but are a finite fraction of the 
mass of the static solution. In analogy to statistical physics this kind of critical 
behavior, where the mass as a function of the parameter p is discontinuous, was 
called type /, whereas the phenomena found by Choptuik, where the mass is a 
continuous function of p are called type II. 

We remark that the study of type II critical phenomena is entirely based on 
classical general relativity, ignoring the fact, that for regions with very strong 
curvature, which necessarily occur in type II critical collapse, the classical theory 
should be replaced by a quantum theory. 

This work is part of a project, which originally aimed at investigating critical 
phenomena of the self-gravitating a model in the presence of a cosmological con- 
stant. From the results of previous work [2] concerning static solutions of the 
model on de Sitter background it was reasonable to expect that these solutions 
would persist when gravity is "turned on". In particular one of the static so- 
lutions constructed in ;2] has one unstable mode, which would therefore be a 
candidate for type I critical collapse^. It would then have been possible to in- 

^In the coupled situation the "end states" could be de Sitter space and Schwarzschild-de- 
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vestigate the coexistence of type I and type II critical phenomena on one hand, 
and the dependece of these phenomena on a dimensionless parameter (see below) 
on the other hand. Therefore we started to investigate existence and stability of 
static solutions in the presence of a positive cosmological constant A. However, 
because unexpected new phenomena emerged, our attention focused on critical 
phenomena in the asymptotically flat situation. The work on static solutions 
with A therefore is rather isolated from the rest of this thesis. 

The main part of this thesis concentrates on type II critical phenomena of the 
self-gravitating SU(2) a model in spherical symmetry (without cosmological con- 
stant). We investigate the phenomena both by evolving one parameter families of 
initial data and doing a bisection search to fine-tune the parameter, and by "di- 
rectly" constructing (i.e. by imposing the symmetry on the equations and solving 
the resulting reduced problem) the relevant self-similar solutions and studying 
their stability properties. 

The motivations for choosing the SU(2) a model as the matter model are, that it is 
a very simple model (in spherical symmetry the field equations reduce to a single 
nonlinear wave equation) and that the theory contains a dimensionless parameter, 
the coupling constant rj. Using dimensional analysis one can expect, that type II 
critical phenomena and the spectrum of self-similar solutions depend strongly on 
this parameter. This expectation is supported by previous work on the limits of 
strong |5l| and weak coupling and [S3], where the solutions at the threshold 
are DSS in the limit rj oo, and CSS in the limit of vanishing coupling. This 
model therefore gives the chance to find out more about the mechanisms that 
are responsible e.g. for the realization of continuous self-similarity as opposed 
to discrete self-similarity. In particular the expected transition of the critical 
solution from CSS to DSS as the coupling is increased is of major interest. 

This work is organized as follows: in Chapter|21the self-gravitating SU(2) a model 
is introduced. We give the basic definitions and equations, that are necessary for 
the work on static solutions, self-similar solutions and critical collapse. The Ein- 
stein equations (with and without cosmological constant) and field equations are 
given with respect to coordinates, that are adapted to spherical symmetry. In 
particular the time evolution code DICE (see App. is based on a character- 
istic formulation of the initial value problem. The coordinates adjusted to this 
formulation (Bondi coordinates) are discussed in this chapter. 

Chapter |31 deals with static solutions of the model in the presence of a positive 
cosmological constant A. We discuss the static equations, investigate the possible 
global structures of solutions and describe the numerical results. The limit rj — > 
Tjmaxi the maximal value of the coupling constant for which solutions exist, is 
carried out with care. 

Chapter m deals with self-similar solutions of the model. We introduce the concept 
Sitter (Kottler) space. 
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of self- similarity and its manifestation in Bondi coordinates. The main part of this 
chapter is dedicated to the numerical construction and linear stability analysis 
of CSS solutions and a discretely self-similar solution. The equations are given 
in adapted coordinates, the numerical methods for constructing the solutions are 
explained, results are discussed and stability of the solutions is investigated. In 
the last section we summarize which self-similar solutions might be relevant for 
the dynamics, either as "end states" or as critical solutions. Furthermore we 
investigate the relation between the DSS solution and the first CSS excitation. 
We put forward the hypothesis that the DSS solution bifurcates from the CSS 
solution in a homoclinic loop bifurcation at some critical value of the coupling 
constant rjc — 0.17. 

In Chapter we describe our results concerning type II critical phenomena. For 
very small couplings the stable CSS ground state gives rise to the formation of 
naked singularities for "intermediately strong" initial data. For these couplings 
we therefore investigate critical phenomena between dispersion and the formation 
of naked singularities. For larger couplings the end states are flat space on one 
hand and black holes on the other hand. As expected, the critical behavior 
strongly depends on the coupling. For large couplings the critical solution is 
DSS, whereas for small couplings it is CSS. At intermediate couplings - in the 
"transition regime" - we find a new behavior, which we call "epsiodic CSS" . This 
behavior is consistent with the hypothesis of the homoclinc loop bifurcation, 
mentioned above. 

Appendix El explains the shooting and matching method, which is used for the 
numerical construction of both static solutions and self-similar solutions. 

In Appendix El we explain the method of discrete Fourier transform, which is 
used for the numerical construction of the DSS solution. 

Appendix O describes the DICE code, that evolves the self-gravitating SU(2) a 
model in spherical symmetry. This code is used to investigate critical phenomena. 

Our conventions concerning curvature quantities and the signature of the metric 
are those used in [TE] (for details see App. |D)). Throughout this work the speed 
of light is set to unity, c = 1. The index notation should be self-explaining. 
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Chapter 2 



The SU(2) (J Model in Spherical 
Symmetry 

In this chapter we give the basic definitions and equations, that will be used for 
the work on static solutions as well as on critical collapse. We give the definition 
of the self-gravitating SU(2) a model (with and without positive cosmological 
constant). On spacetime we introduce two coordinate systems, that are adapted 
to spherical symmetry. For the work on static solutions we choose coordinates 
such that the hypersurfaces of constant time are spacelike, whereas for the work 
on critical collapse Bondi-like coordinates are used. Regularity of the metric in 
these coordinates near the center of spherical symmetry is discussed. Spherical 
symmetry is also imposed on the matter field via the so called hedgehog ansatz. 
We write down the Einstein equations and the field equation in these coordinates 
and give the formulae for the Misner-Sharp mass function as well as the Bondi 
mass and news function for asymptotically fiat configurations. 

2.1 The SU(2) a Model 

The SU(2) 0" model was first introduced into physics by Gell-Mann and Levy j2H] 
in order to describe the meson fields 7r+,7r„,7ro and cr, subject to the condition 
TT^ + TT^ + tTq + cr^ = 1. Geometrically it is an harmonic map^ from spacetime 
(M,5f) to the target manifold {SU{2) ~ S^,G), where G is the standard metric 
on S^. Harmonic maps are well known in mathematics (see e.g. the review 
articles by Eells and Lemaire 1^)- Some of their applications in physics 
are described by Misner j23ISni- In particular Misner points out that harmonic 
maps are geometrically natural nonlinear theories - as are gravity and Yang-Mills 

^To be more precise, as the base manifold is spacetime with a metric with Lorentzian sig- 
nature, the map is caUed a wave map. 
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theory - and could be used to model (the much more involved) nonlinearities in 
general relativity. 

Concerning this work our interest in the SU(2) a model is mainly based on the 
fact that the theory contains a dimensionless parameter, while being a very simple 
generalization of the massless Klein-Gordon field. As will be argued below, both 
static and critical solutions are likely to depend on this dimensionless parameter. 
Therefore we are in the position to study (hopefully) a variety of phenomena 
within a comparatively simple model. 

We start by choosing coordinates on spacetime and coordinates X"^ on the 
target manifold, and denote the map by X, so 

{M,g^,) > {S\Gab) (2.1) 

x^' > X^{x^'). 

The three components X^{x^) of the map are scalar fields on spacetime. The 
harmonic map is defined to be the extremum of the action 

Shm = - y y d^x g^^d.X^'d.X^GABiX). (2.2) 

Geometrically the Lagrangian is the pull back of the metric of the target manifold, 
contracted with the inverse metric on spacetime. fl is the couphng constant of 
the fields. 

Variation of the action with respect to the fields X^ yields the field equations 

UgX^ + g^'' ficiX'') = 0, (2.3) 

where is the wave operator on spacetime = g^^V and T'^q denote the 
Christoffel symbols of the target manifold. For fixed base space (M,5f) Eq. ()2.3|l 
is a coupled quasi-linear system of wave equations for the X^. 

The simplest such system would be obtained, if the target manifold were just 
R, for which the field equation for the single field X would just reduce to the 
Klein-Gordon equation. On the other hand assuming the base manifold to be 
one-dimensional and the target manifold being arbitrary, the system ()2.3p would 
describe a geodesic on the target manifold. In this sense harmonic maps are 
simple geometric generalizations of both the Klein-Gordon field and geodesies. 

In order to incorporate gravity, the action ()2.2j) has to be supplemented by the 
Einstein-Hilbert action 

Stot = - 2A) + Shm, (2.4) 

where TZ denotes the Ricci scalar of spacetime, k contains the gravitational con- 
stant G hj K = SttG. In addition we have introduced a cosmological constant A, 
which we will always consider positive in this work. 
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Variation of the total action ()2.4|1 with respect to the spacetime metric g^j^u yields 
the Einstein equations 

G^iv + g^iA = nT^iy, (2.5) 

where T^,y is the stress energy tensor of the a model, obtained from varying the 
matter part of the action ()2.2j) with respect to the spacetime metric 

T,u = fl (^V^X^V.X^ - ^(^M.IV.X^V'^X^)^ Gab{X). (2.6) 

This stress energy tensor satisfies the weak, strong and dominant energy condi- 
tions, as defined in |39j. The weak energy condition (WEC) requires T^^w'^w^ > 
for all timelike w corresponding to a positive energy density for all observers. The 
strong energy condition (SEC) results from the "timelike convergence condition" 
Rfj_yW^w'^ > for all timelike w, which for A = translates into a condition for 
the stress energy tensor T^yW^w^ > (l/2)w^w^T. Finally the dominant energy 
condition (DEC) consists of the weak energy condition plus the requirement that 
T^^^Wy is non-spacelike for all timelike w. An equivalent condition is, that the 
components in any orthonormal basis satisfy > \T^'^\ for all /, J. 

For the WEC we choose an orthonormal basis {e7-}/=o,i,2,3 with the coordinate 
representation ej = e^d^. We consider the expression T^yC^eQ = Tgg, which can 
be written as 

T>.Aeo = l{X,eo)^{X,eofGAB + ^{X,ei)^{X,e,f6'^GAB > 0, (2.7) 

as Gab is Riemannian and therefore the above expression is a sum of positive 
(or vanishing) terms. By Xj,v we denoted the push forward of the vector v 
from spacetime to the target manifold. As the above relation is valid for any 
orthonormal basis, the relation T^yW^w'^ > is satisfied for all timelike w. 

For the SEC again we choose an orthonormal basis as above. Inserting this into 
the expression T^yCQC^ + (1/2)T one finds that 

Again this is valid for all orthonormal bases, and therefore the matter field sat- 
isfies the SEC borderline, as the massless Klein-Gordon field does. 

Finally to check the DEC we work with the components of T^^, in an orthonormal 
base. We have 

T°° = ^^ = ]^{X,e^)\X,e,)^GAB + \{X,e,)\X,e,)^5'^GAB. 
= -Tq. = -(Xeo)^(X,e,)^G'^B 

= J]. = {X,e,)^{X,ejfGAB - ^S,, [-(X,eo)^(X,eo)^ + 5''(Xe,)^(X,Q)^] G 
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As Gab is positive definite, we have \2G{v,w)\ < G{v,v) + G{w,w). Therefore 
|T*| < T^*^. One also easily sees, that |T**| < (without summation over i) 
and that \T'^^ < for i ^ j. Again the tetrad was chosen arbitrarily, so the 
relations are valid in any orthonormal base. 

Note that the coupling constants fl and G only enter the equations ()2.5|1 as the 
product 7] := In units where the speed of light is set to unity c = 1, G 

has dimension of length/mass and has dimension of mass /length, so their 
product is dimensionless. Therefore the only scale in the theory is tied to the 
cosmological constant A, which has dimension of 1/ length?. As A merely sets 
the scale of the theory, its actual numerical value can be eliminated from the 
equations and only it's sign matters. 

For fixed A on the other hand the presence of the dimensionless product of cou- 
pling constants rj := provides a one parameter family of physically different 
theories. 

For A = 0, the theory is scale invariant. This implies immediately (as will be 
explained in Sec. 12.1.11) that for A = this model does not admit asymptotically 
fiat soliton solutions 6^ . On the other hand, as described in Sec. I4.1.H| the scale 
invariance for A = is necessary for the existence of self-similar solutions. 



2.1.1 Non-existence of Asymptotically Flat Soliton Solu- 
tions and Static Black Holes for A = 

Using a scaling argument, it is easy to see that the system ()2.Hj) and ()2.5|1 does 
not admit static, globally regular, asymptotically fiat solutions (so-called soliton 
solutions) apart from the trivial solution, which is Minkowski spacetime with 
vanishing matter field. Assume there exists a static solution {g^y,X^) to the 
system fl2.H|l and ()2.5|1 . We denote the timelike Killing vector by ^ = dt and choose 
coordinates in the hypersurfaces S orthogonal to ^. For regular solutions S is 
topologically R^. Denoting the induced metric by hij we can write ds"^ = —Ndt^ + 
hijdx^dxK Staticity in these coordinates manifests itself in the independence of 
the metric functions N and hij and the fields X"^ of the time coordinate t^. 

For static solutions, due to the existence of the Killing vector C,, one can define 
the energy of the field 

E = - J Vh d^x n^ C = - J d^x N f °. (2.10) 
s s 

Asymptotic fiatness guarantees the existence of this integral. E is a conserved 
quantity, i. e. independent of the hypersurface S, which follows from V^(^'^Tj^) = 

^For implications of a spacetime symmetry on the fields see Sec. 12.2.31 
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0. Furthermore E is positive definite as the lapse is positive, and the energy 
density — Tq is positive as well. Moreover Tq vanishes iff (X^,ej) = (or equiva- 
lently X is constant), i.e. if the map is trivial. 

Since the model for A = is scale invariant, there exists a one-parameter fam- 
ily of solutions (^{gfj,u)x{x'^) '■= g^lu{^x^), ^x{x^) '■= X^(Aa;*)j. The corresponding 
energies Ex are obtained from E hj Ex = As a static solution extremizes 

the energy ()2.1()j) we must have dEx/d\\x=i = 0, so E has to be zero and therefore 
the matter field has to vanish, = const. We are left with a regular, static, 
asymptotically flat solution to the vacuum Einstein equations, which has to be 
Minkowski spacetime due to the theorems by Lichnerowicz [H^. 

It was also shown 02] ? that this model with A = neither admits non-trivial 
static asymptotically fiat black hole solutions. 

So if one is interested in static solutions of this model, the cosmological constant 
is essential, or as we shall see in Sec. IH.2.2| the assumption of asymptotic flatness 
and S being has to be dropped. 



2.2 Spherical Symmetry 

An isometry is a diffeomorphism $ from (M, g) to (M, g) ,which maps the metric 
to itself, i.e. 

^*g = g. (2.11) 

For a one parameter family of such diffeomorphisms $a, with $o = "id one can 
define the generator ^ = d^x/dM>y=o- Eq. ()2.1H) then can be formulated in terms 
of the Lie derivative 

C^g = 0, (2.12) 

or equivalently 

^(m;.) = 0. (2.13) 

Such a family of isometrics leaves the curvature tensors invariant (this can be 
seen in taking the analogous steps as in Sec. 14.11) . in particular 

(AG),. = 0. (2.14) 

A spacetime is said to be spherically symmetric, if it admits the group SO (3) as a 
group of isometrics, acting on spacelike two-dimensional surfaces (See e.g. |H^). 
The group acts transitively but not simple transitively: as the group is three 
dimensional and the spacelike surfaces are only two-dimensional, it has SO (2) 
as an isotropy group. The orbits of the group are surfaces of constant positive 
curvature. 
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It can be shown that the metric of a spherically symmetric spacetime can be 
written as the warped product 



ds^ = dr^ + R\t') {d9^ + sin2 9 d^^) , (2.15) 

where the coordinates 6 and ip are coordinates on - the orbits of SO (3) - with 
the usual range < < tt and < < 27r. rfr^ denotes the line element of a 
two-dimensional Lorentzian manifold (with coordinates r*), and -R(t*) is related 
to the area of the orbits of S0(3) hj A = AirR"^. 

For this work we use two different choices of coordinates for the Lorentzian two- 
surfaces: for the work on static solutions (in the presence of a positive cosmolog- 
ical constant), we choose orthogonal coordinates (t,p) 

ds'^ = -A{t, p)dt^ + B{t, p)dp^ + R^{t, p)dn'^, (2.16) 

where dQ^ = {dO"^ + sin^ 6' dip^). Usually these coordinates are further restricted 
by choosing the area of the orbits of SO (3) as a coordinate. But this is only 
possible in regions, where V^i? 7^ 0. As some of our static solutions will contain 
hypersurfaces where the gradient V^-R vanishes, we choose a different gauge: we 
set l/B = A=: Q, so 

ds" = -Qit,p) de + + R\t,p) dn\ (2.17) 

For the investigations on critical collapse, we work with Bondi-like coordinates, 
i.e. we foliate spacetime by outgoing null cones u = const, which emanate from 
the center of spherical symmetry {R = 0) and parameterize these with the area 
of the two-spheres r := A/4tt, so we get 

ds^ = -e^i^M du (^{u, r) du + 2dr^ + r^dn\ (2.18) 

The normal vectors to the hypersurfaces u = const are null, as — V^m = (— 1, 0, 0, 0) 
and — V^M = (0, —g^^, 0, 0), so V^mV^m = 0. Furthermore, they generate affinely 
parametrized geodesies, as 

V'mV^V'm = V'mV^V^m = ^V^ (VauVu) = 0. (2.19) 

The areal coordinate r parameterizes these null geodesies, but is not the affine 
parameter. 

Ingoing radial null geodesies are given as the solutions of 

drju) ^ V{u,r{u)) 
du 2r(u) ■ ^ ■ ^ 
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We normalize u to be proper time at the origin. We require j3{u,r = 0) = 0. 



Regularity at the origin then enforces 



M, r 



0) = 1 (see Sec. m^ . The 



connection between these "Bondi-like" coordinates and Bondi coordinates for 
which PB{uB,r = oo) = is given by a coordinate transformation u Ub{u), 
with 



in particular at infinity we have 
duB 



du 



Therefore 



PB{uB,r) 
— I iuB,r) 



with H(u) 



oo 



(3{u,r)-H{u) 



r 



B 



V 



— I (M,r)e 
r 



-2H(u) 



(2.21) 



(2.22) 



(2.23) 



For further use we give the square root of the determinant of the metric 
e^^r^ sin 9 and the inverse metric 



( ° 



e 








2PV Q 














'-g 



(2.24) 



2.2.1 A Regular Center of Spherical Symmetry 

Clearly the coordinates ()2.17|) and ()2.18|) break down at a center of spherical 
symmetry = 0, r = respectively. Apart from the vanishing volume of the 
two-spheres, which is well known from polar coordinates in flat space, the metric 
functions Q, R and /3, ^ have to satisfy additional regularity requirements, if 
one asks for a regular center of spherical symmetry. By definition the metric is 
regular (C*^), if its application to regular {C^) vector fields yields regular iC^) 
functions on the manifold. The easiest way to examine this is to switch to regular 
coordinates close to the center. 

We start with the coordinates ()2.17j) . First we fix the origin of the coordinate p to 
be at the center of spherical symmetry, R{t,p = 0) = 0. We choose coordinates 
{t, x,y,z), connected to (t, p,9,Lp) by x = p sin 6^ cos (p, y = p sin 9 simp, z = p cos 9 
and t = t. We assume (t, x,y, z) to be regular coordinates in the vicinity of 
p = 0. A function then is regular, if it can be written as a regular function of 
these coordinates. By specifying the above coordinate transformation we also 
have implicitly assumed that the coordinate function p has special regularity 
properties, namely p itself is not regular, but any even power thereof is. First 
note that g{dt, dt) = —Q, which is regular if Q is a regular function of (t, x, y, z). 
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In other words, Q has to be a regular function of t and p^. Second we consider 
the sum of the spatial components of the metric with respect to the regular 
coordinates, 

g{d,, d,) + g{dy, dy) + g{d,, d,) = ^ + 2^. (2.25) 

This shows, that R/p has to be a regular function of t and p^. Having again a 
look at g{dz,dz) , 

1 sin^^ 

g{d,, d,) = cos' e- + —^R\ (2.26) 

Q p 

we see, that the only possibility for this expression to have a regular limit p — > 0, 
is 

lim^ = l, (2.27) 

P^O p2 

SO Q{t,Q) = 1/R'(t,0), which we can choose without loss of generality to be 
1. Therefore near a regular center of symmetry the metric functions behave as 
follows 

R{t,p) = p + 0(p3) 

Q{t,p) = l + 0(p2). (2.28) 

For the Bondi-like coordinates ()2.18|1 we define in analogy the coordinates t = 
u+r, X = rsinO cos (p, y = r sm9 sin <yj and z = r cos 9. We have g{dt, dt) = —e^^^y 
and 

g{d., d.) + g{dy, dy) + g{d,, d,) = 2- e'^"^ + 2e'^. (2.29) 

From this it follows that the metric functions /? and — have to be regular functions 
of t and r^. Looking at g{dz, dz) 

g(dz, dz) = e'\2 - ^) cos^ e + sin^ ^, (2.30) 

we see that lim e^^(— — 2) = 1 is a necessary condition for regularity. P{u,r = 0) 

has already been chosen to be unity, so ^(t, r = 0) = 1. It remains to transform 
these functions of t and back to functions of m = t — r and r. The first terms 
in a Taylor series expansion give 

f3{u,r) = 0{r') 

j{u,r) = l + 0{r'). (2.31) 

Note however, that /3 and ^ if expanded in u — uq and r don't solely contain even 
powers of r. The first non- vanishing term with an odd power in r is e.g. the term 
l3"iuo,0)r\ 
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2.2.2 The Mass Function 



In spherical symmetry one can define the Misner-Sharp mass function jSZ| m{T^), 

by 

2?77/ 

l- — = V^R V^i?, (2.32) 

rC 

(for a recent description of the properties of this function see the article by Hay- 
ward IHj) which gives 

m{t,p) = f (l + R'{t,p) - Q{t,p){R'{t,p)r'^ (2.33) 

for the coordinates ()2.17p . Its interpretation for static solutions in the presence 
of a positive cosmological constant will be described in Sec. 13.1.41 

For the Bondi-like coordinates ()2.18|) we have 

m{u, r) = ^ (^1 - e-2^("''') j {u, r)^ . (2.34) 

For an asymptotically fiat spacetime, the Bondi mass is obtained by taking the 
limit r — > oo of m{u, r) along the null hypersurfaces u = const, so 

"^Bondi (m) = lim m(M,r). (2.35) 
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In general, due to radiation (of the matter fields only in spherical symmetry), the 
Bondi mass will decrease with retarded time u. We give an explicit formula for 
the mass loss in Sec. 12.3.11 

The formation of an apparent horizon is signalled by the vanishing of the expan- 
sion of outgoing null geodesies 0+ = 0. 0+ can be expressed as the Lie derivative 
of the area A of 2-spheres with respect to the tangent to outgoing null geodesies 
l'^, divided by the area: 0+ = {Ci_^_A)/A. 

For the Bondi-like coordinates, we have already seen that — V^m is tangent to 
the outgoing radial null geodesies {u = const). The area of the 2-spheres is given 
by y4 = 47rr^, so we have 

0+ = = (2.36) 

An apparent horizon in these coordinates therefore manifests itself by /5 — *• oo. 
The breakdown of these coordinates at an apparent horizon is due to r (in con- 
nection with u) ceasing to be a good coordinate. As stated above the areal 
coordinate r parameterizes the null cones u = const. This is possible as long 
as V„r 7^ 0, which evaluates precisely to g'^^ ^ 0, which is then violated at the 
apparent horizon. 
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2.2.3 The SU(2) a Model in Spherical Symmetry 

Imposing a symmetry on spacetime also requires some symmetry properties for 
the matter. The stress energy tensor T^y has to be invariant under the isometry 

C{r^"' = 0, (2.37) 

otherwise the symmetry of spacetime would be incompatible with the Einstein 
equations. 

For the SU(2) a model we write the stress energy tensor ()2.fi|l as 

T^, = iX*G)^, - ^g,,g''^iX*G)^r. (2.38) 
The Lie derivative of this expression then gives 

£^T^, = C^{X*G)^, - ^g^,g''^C^{X*G)^r, (2.39) 
if ^ is a Killing vector field. 

As the Lie derivative commutes with the pull back 

C^iX*G)^, = iX*iCx,^G))^,, (2.40) 

the requirement ()2.37p is satisfied if 

Cx,(Gab = 0. (2.41) 

This means that either X^^ = or the Killing vector field ^ is mapped to a 
Killing vector field on the target manifold. 

Applying this to spherical symmetry there are essentially two possibilities to make 
the map spherically symmetric. First, if none of the fields X^ depends on the 
angular variables 6 and if then the Killing vector fields C,i of spherical symmetry 
on spacetime would be mapped to the zero vector field at the target manifold. 
This way one would deal with three fields X^{t^). 

The second possibility, which is chosen in this work, uses the symmetry of the 
target manifold: as the metric Gab is the metric of constant curvature on 
it also admits SO (3) acting on 2-spheres as a group of isometrics. We choose 
coordinates (0, B, $) such that the line element is given by 

dsls = d(f)^ + sin^ (f){dQ^ + sin^ 9 d<^>^). (2.42) 

Obviously the vector fields 

51 = sin $c}e + cot 6 cos $9$, 

52 = — cos $(9e + cot 6 sin $5$, 

53 = -d^ (2.43) 
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are Killing vector fields on the target manifold. We demand now, that the corre- 
sponding Killing vector fields on spacetime are mapped to their counterparts 
on the target manifold, i.e. 

= S,. (2.44) 
This is achieved by the so called hedgehog ansatz 

(j){x^') = 0(r^), eix") = e, ^x^') = 0. (2.45) 

This way the field equations ()2.3|) decouple into a nonlinear wave equation for 
0(r*) and two equations for G and $, which are satisfied identically. 

In order to examine regularity of the map at the center of spherical symmetry, we 
again work with Cartesian coordinates X = (p sin B sin $ = sin 6 cos (f = (0/r)x, 
Y = 0sinBcos$ = (psinOsimp = {(j)/r)y, Z = 0cos0 = (pcosO = {(j)/r)z 
(for the coordinates ()2.17|) r is replaced by p). The fields {X,Y,Z) are regular 
functions on spacetime, iff 0/r (0/p respectively) is a regular function. Therefore 
close to the origin we get the expansions 

0(t,p) = p(l + 0(p2)) (2.46) 
(f){u,r) = r(l + 0(r^)). (2.47) 

for fixed time t or fixed retarded time u. Again for fixed t, is a regular function 
of whereas for fixed u odd powers of r appear. 

In particular this means that for all times t (u) the origin is mapped to a single 
point on S'^, which is the north pole as defined by ()2.42|) . 



2.3 The Einstein Equations and the Field Equa- 
tion 

With the hedgehog ansatz ()2.45|) the field equations ()2.3|) reduce to the single 
wave equation 

^ , sin(20) ,sin(20) , / . x 

□.0 = -^ (^^resp.), (2.48) 



where 0^0 reads 



□^0 = ^(^-dt{^d,) + d,{R'Qd,)]<p (2.49) 



Dg<P = e-2'^(^(^^+(^^^^9,-^5„-2a„a, + ^9,,^0, (2.50) 
for the coordinates ()2.17|) . ()2.18|) respectively. 
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For the work on static solutions (see Chapter ^ the following combinations of 
the Einstein equations + A6^ = t^Tj^ will turn out to be convenient: the 
combinations (*) + (p — 2(g), 



2g(p,t)- 



Qip.t) , -2-'-§^ + 2Q{p,t)R'{p,tr 



-2 sin((/)(r,t))2 0(r,t)2 



Q"{p,t) 



2r] 



R{r,ty 



Q{r,t) 



+ Q(r,t)0'(r,t)2 , (2.51) 



2 {R{p,t) + Qip,tfR"{p,t)] / w ^2 

-2v t^+QiP,m'ip,tr 



Q{p,t) R{p,t) 



Qip,t) 



(2.52) 



of the Hamiltonian constraint {j) and the time evolution equations (p and (g), 



the time evolution equation 



A - R{p,t)-' + 



R{p,tf 



2R{p,t) 



Q{p, t) R{p,t) 

Q(p, t)2 R{p, t) Q{p, t) Rip, ty Q{p, t) R{p, t) 



+ 



^ Q'ip,t)R'ip,t) ^ Qip,t)R'ip,tf 



Rip,t) Rip,ty 

-2sin(0(p,t))2 , 0(p,t)2 



R{p,ty 



+ 



Q{p,t) 



+ Q{p,t) (j)'{p,tf 



(2.53) 



and the momentum constraint (* ) 



-R{p, t) Q'ip, t) ^ Q{p,t)R'{p,t) ^ 2 R'ip, t) 



Q(p, t)mip, t) Qip, t)mip, t) Q{p, t)R{p, t) 



where ' = dp and ' = dt- Of course = Gg and all other components vanish 
identically. (See ChapterElon the structure of these equations in the static case.) 

For the work on critical collapse in the coordinate frame ()2.18|1 the nontrivial 
Einstein equations split up into the hypersurface equations (the {rr} and {ur} — 
{V/2r){rr} components of Gf^u = i^Tfj^u) 



v = 

and the subsidiary equations —E^ 
and Eat, = Gm - SnToa: 



e2^(l-2r/sin(0)2), 



(2.55) 
(2.56) 



E,, 



2V/3 -V + 27]/' 



r 



0, 



(2.57) 
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V{r(5" - 13') + rl3'V' + -rV" - 2r^$' 
-r/ rVf--0' + +20) =0, 



(2.58) 



r 



where ' = dr and ' = du- 

The contracted Bianchi identity V^jG^^ = together with V ^Tj^ = for solutions 
of show, that the system (jOHI), and ^IT^i is sufficient to evolve 

the Einstein a model: The {66} component of the Einstein equations, Eqq = is 
satisfied, if Eqs. (jTiH|l . (P3?)|l and (E37j) are satisfied. This follows from 

the "r-component" of the Bianchi identity: Ege can be expressed as an algebraic 
combination of the other components of the Einstein equations and derivatives 
thereof. 

Furthermore the "w-component" of the Bianchi identity reads (r'^Guu)' = {f^VGur)'— 



e^^r'^du{e-^^Gur) + Ire^^GrrduiVe-'^^). Assuming, that Eqs. (EHI), and 
fl2.56p are satisfied, then {r'^{Guu — f^Tuu)') = 0. Therefore r'^{Guu — i^Tuu) = f{u)- 



So if Guu = i^Tuu at some hypersurface r = const then this equation is satisfied 
everywhere. Now the regularity conditions at the origin ensure that this equation 
is satisfied at the origin r = 0, so the {uu} component of the Einstein equations 
is satisfied everywhere if the field equation ()2.48|) and the hypersurface equations 
()2.55p and ()2.56|) are satisfied. 

So the characteristic initial value problem we deal with consists of the system 
()2.48|) . ()2.55|) and ()2.56|) together with the initial conditions (f){u = 0,r) = 4>{r). 
The numerical treatment of this characteristic initial value problem is described 
in Appendix O 

2.3.1 Mass Function, Bondi Mass and News Function 

We close this section by giving explicit formulae for the mass function m(M, r) 
(I2.34|) and the mass loss at infinity rriBondiiu) in terms of the matter field. 

The mass function m(u, r) can be rewritten as an integral over outgoing null rays 



r 



by the trivial identity m{u^r) = J m'{u,r)dr. m'{u,r) is given by 






(2.59) 



Using the hypersurface equations ()2.55j) and ()2.5(ij) we obtain 




(2.60) 



and therefore 



r 




(2.61) 
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This expression together with ()2.34j) will serve as an accuracy test for the numer- 
ical code (see Appendix [n)). 

In order to derive a formula for the mass loss at infinity in terms of the matter 
field we have to look at the behavior of the metric functions and the field at 
infinity. The Bondi mass ()2.35|) is finite if (3 = H{u) + ai{u)/r + 0(l/r^) and 
Y. = (,2H{u) + 0(l/r2)). From (piH]) follows further, that 0(M,r) = 

+ 0(l/r^). Using the hypersurface equation ()2.55|) we find that ai{u) = 0, 

so 

(3{u,r)=H{u) + 0{l/r'^). (2.62) 

Inserting these expansions into the formula for the Bondi mass fj2.35p we get 

hi{u) = -2m£,o„die^^("\ so 



-{u, r) = e2^(") ( 1 - + 0{l/r'). (2.63) 



The derivative of the mass function m{u, r) with respect to retarded time u is 
given by 

m{u, r) = -r^r^e"^^ ( 0^ (j)'^] , (2.64) 



r 

where one of the subsidiary Einstein equations ()2.57|) has been used. This gives 
for the mass loss at infinity 

rfiBondiiu) = -r] c?(u)e-2^("). (2.65) 

Clearly ntBondi < corresponding to the energy, that is radiated away to infinity 
and therefore lost. The expression 

N{u) := c?(u)e-2^(") (2.66) 

is called the news function. 
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Chapter 3 

Static Solutions of the 
Self-gravitating a Model in the 
Presence of a Cosmological 
Constant 



3.1 The Static Equations 

In this chapter we investigate the question of existence and stabihty of static so- 
lutions to the self-gravitating SU(2) a model with positive cosmological constant 
in spherical symmetry. This work is motivated by previous work |2] on static 
solutions of the model on fixed de Sitter background and by the work of Volkov 
et al. [ZSj and Brodbeck et al. [15j, who considered the Einstein Yang- Mills system 
with positive cosmological constant. 

The choice of coordinates as well as the introduction of a gauge invariant quantity 
for the stability analysis follows [T3] . In addition we put some emphasis on 
examining the role of the Killing horizon. Our results closely resemble those of 
ff5| IT5j with the only difference, that in the limit of maximal coupling (see Sec. 
IH2.2|) the system is scale invariant. The results on static solutions are summarized 
in IHl. 



3.1.1 Staticity 

A spacetime (M, g) is called stationary, if it admits a timelike Killing vector field 
^. If this Killing vector field is in addition hypersurface orthogonal, then the 
spacetime is called static. Hypersurface orthogonality is given if the Frobenius 
condition 

^[aV^e.] = 0. (3.1) 
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is satisfied. 



Consider now a static, spherically symmetric spacetime (see e.g. (ZHj)- If the 
static Killing vector field ^ is unique, it has to be invariant under the action of 
SO (3) (As the composition of two isometrics is again an isometry, an element of 
SO (3) maps ^ to a Killing vector field. Furthermore the norm is left invariant, 
so ^ is mapped to a timelike Killing vector field. If ^ is unique, then SO (3) 
has to map C, to itself). From this it follows, that ^ must not have components 
tangential to the orbits of S0(3), as the only vector field on S^, which is left 
invariant under the action of SO (3), is the zero vector field. Therefore ^ can be 
written as ,^ = ^^{t,p)dt + ^''(t, p)dp. We are still free to choose the coordinates 
(t, p) such that ^ = dt- With this choice the metric ()2.17|) in the presence of a 
hypersurface orthogonal Killing vector field ^ = dt reads 

ds^ = -Q{p)df + + R\p)dn\ (3.2) 

Q{.p) 

Stationarity of spacetime again extends to the field via the Einstein equations. 
In order to satisfy ()2.37|) we set = 0(p). 



3.1.2 The Static Equations 

Setting all time derivatives in Eqs. ()2.48p and ()2.5H) - ()2.54|) to zero, we find that 
the momentum constraint ()2.54|) is satisfied identically. The field equation ()2.48|) 
and the combinations 1)2.511) and ()2.52|1 yield the following autonomous system 
of coupled, nonlinear, second order ODEs 

{QR^cj)')' = sin 20, (3.3) 
{R^Q'y = ~2kR\ (3.4) 
R" = -r/i?0'2. (3.5) 

Furthermore Eq. ()2.53p (multiplied by R^) is a first integral of the above system: 
2r/ sin^ + R^{k - r/ Q0'2) + RQ' R! + Qi?'^ -1 = 0. (3.6) 

From the p component of the contracted Bianchi identities we have {y^R'^G'^' = 
i?^(x/Q)'G* + y/Q{R^yGl). Assuming Eqs. (jSISl)-(|S31) to be satisfied, we have 
QR?(G''p — KTp = const. So if ()3.6|) is satisfied at some hypersurface p = const it 
is satisfied everywhere. Now the conditions ()2.28j) and ()2.4(j)l for a regular center 
of spherical symmetry yield, that Eq. ()3.(i|l is satisfied at p = and therefore it 
is satisfied everywhere. 

As the cosmological constant A - if non-zero - sets the length scale of the theory 
it can be eliminated by switching to the dimensionless variables p = a/Ap and 
R = ^/AR. Or in other words: setting A to unity (say) in the above equations 
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means that all quantities which have dimension of length are measured in units 
of VA. 

We are interested in solutions of ()3.3p - (j3.5p which have a regular center of spher- 
ical symmetry = (at p = 0) and a Killing horizon Q = at some finite 
distance from the origin. (See Sec. 13.1.41 for the global structure of such space- 
times). The reason why we look for solutions with a horizon is the following: 
turning off gravity [r] = 0) Eqs. ()3.3j) - ()3.5|1 (in combination with the regularity 
conditions at the origin) describe the de Sitter spacetime (see Sec. I3.1.3|) . which 
has a cosmological horizon, and we don't expect the global structure to change 
when gravity is turned on slightly [rj small and positive). Furthermore we can 
rule out the following possibilities: integrating outward from a regular center 

1. the static region "ends" in a singularity, 

2. the static region has a second (regular) pole R = 0^, 

3. the static region persists up to spatial infinity. 

In principle the first case could be produced easily by choosing an arbitrary value 
(f)'{p = 0) at the regular center (See Sec. 13. 2^ . Nevertheless this case is of no 
interest here and is therefore discarded. Avoiding this scenario means, that we 
have to set up a boundary value problem enforcing one of the other cases. 

The second case is impossible for A > 0: rewriting equation ()3.4|) as an integral 
equation ()3.19j) . one sees, that Q' diverges, whenever R goes to zero for a second 
time. 

That the last case is impossible, can be seen as follows: assume the static region 
extends to infinite values of p, i.e. Q{p) > for all p > 0. Assume further, that 
R'{p) > for all p > (this second assumption is necessary, as R'{p) is strictly 
monotonically decreasing (see Eq. ()3.20|) ). so once it gets zero R decreases until it 
eventually becomes zero, which in turn would lead to the cases 1 or 2). Then one 
can show that for p large enough Q' is bounded from above by Q'{p) < —const/ p. 
This means that Q has to cross zero at some finite value of p, so case 3 is ruled 
out as well. 

We will discuss the implications of Eqs. ()3.3|) - ()3.6p for A = in Sec. ()3.2.2|) . 
3.1.3 Exact Solutions 

In this section we give some simple exact solutions of Eqs. ()3.3|) - ()3.5|) . which will 
be important for the full spectrum of solutions (See Sec. I3.2|l . as they arise in 
certain limits. 

^Note that this situation is not ruled out if A = 
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We start with the solution obtained by setting = 0. Then R"{p) = and 
therefore (together with the gauge choices -R(O) = 0, R'{0) = 1) we have R{p) = p- 
Eq. (1131) then gives Q'{p) = — 2Ap/3 and Q{p) = 1 — Ap^/3. Together we have 

0(p)=O, Rip)=p, Q(p) = l-^, (3.7) 
which is the de Sitter spacetime in the static frame. 

Another simple solution is obtained by setting = 7r/2. This violates the reg- 
ularity condition for the field at the origin. Again we have R{p) = p, but Q{p) 
cannot fulfill the regularity requirement Q{0) = 1/R'{0) at the origin any more, 
instead Q{p) = 1 — 2// — Ap^/3, so 

0(p) = ^, i?(p)=p, Qip) = l-2r]-^. (3.8) 

This solution behaves like the de Sitter spacetime for large values of p, but has 
a conical singularity at the origin. In the limit of vanishing coupling 77 = 0, 
where spacetime is de Sitter, the solution = 7r/2 has diverging energy density 
at the origin, but finite total energy (as measured between origin and horizon) 
maximizing the energy of all regular static solutions that exist on fixed de Sitter 
background and can be viewed as a "high excitation" limit of this spectrum (see 
0). 

For A = we obtain a solution, which will be of interest in the limit of maxi- 
mal coupling 1] — i> rjmax (see Sec. I3.2.2|l . Eq. ()3.4|1 together with the regularity 
conditions at the origin give Q{p) = 1. The remaining equations can be solved 
analytically for rj = 1 to give the static Einstein universe: 

0(p)=p, i?(p) = sinp, Q(p) = l, r/ = l. (3.9) 



3.1.4 Horizons and Global Structure 

Clearly the coordinates ()3.2|) break down at the horizon Q{ph) = 0. In order 
to examine the global structure of spacetime we temporarily switch to outgoing 
Eddington-Finkelstein coordinates {u, r, 6, 0), where the retarded time u is given 
by 

(3.10) 

The metric fj3.2|) then reads 

ds"^ = -Q{p)du^ - 2du dp + R{pfdVt^. (3.11) 

Note that the coordinates ()3.1H1 cover only half of the maximally extended space- 
time. In order to cover all of spacetime one would have to switch to "Kruskal-like" 
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double null coordinates, as is done e.g. in [ZSj- In the following, we will simplify 
our discussion by only talking about the Killing horizon contained in the portion 
of spacetime covered here. All statements made can be extended trivially to the 
complete spacetime and in particular the second component of the horizon by 
time reflection. We also remark that all solutions have the topology x R. 

The static Killing vector field is du = dt, The metric fj3.11|) is regular if Q{p) and 
R{p) are regular functions, except when R = 0, which corresponds either to the 
usual coordinate singularity of spherical symmetry, which has been discussed in 
Sec. (12.2.11) or to a spacetime singularity, as discussed below. 

As described in Sec. 13.1.11 we assume the Killing vector field 9„, to be timelike 
in some neighborhood of the center R = 0, i.e. g{du,du) = —Q{p) < 0. From 
the discussion in Sec. 13.1.21 it is clear, that Q{p) has to go to zero at some finite 
distance from the origin, Q{ph) = 0. Furthermore Q{p) changes sign there (as 
Q'{p) < the degenerate case of Q'{ph) = is impossible). In regions, where 
the norm of the Killing vector is positive, i.e. Q{p) < 0, spacetime is dynamic. 
The slices of constant time p are of topology R x S^. They are generated by the 
Killing vector fields du and and are therefore homogeneous. Such regions thus 
correspond to Kantowski-Sachs models. 

The hypersurface p = pu separating the static and dynamic regions, is char- 
acterized by Q{ph) = 0, and is a null hypersurface (as V^mVm = and 
V^mVo-V^m = 0). As the Killing vector field du is null on and tangent to this 
hypersurface it is a Killing horizon. 

In order to characterize this horizon further, we use the concept of trapping hori- 
zons introduced by Hayward jlO] . In general for asymptotically flat (and asymp- 
totically de Sitter) spacetimes the asymptotic region can be used to classify an 
event horizon as a black hole horizon (or as a cosmological horizon) . Furthermore 
for the (local) concept of an apparent horizon the asymptotic region is needed to 
define inward and outward directions. 

The concept of a trapping horizon is based solely on the local behavior of null 
congruences. As in Sec. 12.2.21 the expansions Q± of out- and ingoing null rays 
emanating from the spheres R = const are defined as 

0± = ■^'^±R'^ (3-12) 
where C± denotes the Lie derivative along the null directions 

l+ = dp and l_ = 2du-Qdp (3.13) 

respectively, so 

0+ = 2^ and 0_ = -2Q^. (3.14) 

R R 

The expansions can be used to define a trapped surface in the sense of Penrose [02] 

as a compact spacelike 2-surface {R = const) for which 0_0_|_ > 0. If one of the 
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expansions vanishes, the surface is called a marginal surface. For a non-trapped 
surface 0_ and 0+ have opposite signs, which is used to have a local notion of 
inward and outward: the direction in which the expansion is positive is called 
outward, and in which the expansion is negative is called inward. 

The closure of a hypersurface, which is fohated by (non-degenerate^ and non- 
bifurcating^) marginal surfaces is called a trapping horizon. At a trapping horizon, 
the behavior of the non-vanishing expansion can be used to distinguish between 
future trapping and past trapping horizons, i.e. supposing 0-1^ = (as is the 
case in our example) then the horizon is 

future trapping if 0+|j? < 0, 
past trapping if 0+\h > 0. (3.15) 

The case Q+\h < means that in the trapped region both outgoing and ingoing 
future directed null rays converge whereas in the case Q+\h > 0, the null rays in 
the trapped region converge if past directed. 

Furthermore the change of the vanishing null expansion 0_ in direction of 
leads to a classification of outer and inner trapping horizons, i.e. the horizon is 

outer trapping if C^Q_\h < 0, 

inner trapping if >C+0_|// > 0. (3.16) 

The meaning of inner and outer are the following: For £_|_0_|// < starting from 
the non-trapped region, where the directions inward and outward are defined as 
above, one has to move inward to approach the horizon, or in other words, the 
horizon is an outer boundary for the trapped region, whereas for >C+0_|ij > 
the horizon is an inner boundary of the trapped region. 

A future outer trapping horizon therefore is suited to describe a black hole hori- 
zon, a past outer trapping horizon would describe a white hole horizon. Inner 
horizons on the other hand describe cosmological horizons, as they occur for 
example in the de Sitter spacetime. 

Our situation is as follows: on the Killing horizon Q = 0, 0_|q=o = while 
Q+Iq^o = 2R'h/Rh ^ and £+0_|q=o = -2Q'hR'hIRh ^ (except when also 
E!h — 0, which we exclude for the moment). 

As i? > and Q' > for p > 0, the character of the KiUing horizon is given by 
the sign of R'^. For i?^ > we have Q+\h > and £+0_|h < 0. The Killing 
horizon therefore is a past inner trapping horizon, and therefore a cosmological 
horizon. For R'^ < on the other hand, wc have 0+|ij < and Cj^Q-\h > 0, 
which corresponds to a future outer trapping horizon, and therefore to a black 
hole horizon. 

Hov e_ = : £+e_ ^ 

Hor e_ = : e+ ^ 
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If i?' = anywhere in spacetime, then both expansions vanish. However, this 
does not mark a trapping horizon, as the product 6+0_ does not change sign (9+ 
as well as 0_ both cross zero). Instead this means, if it occurs in a non-trapped 
region, that the meaning of inward and outward directions are interchanged. 

Using the expansions 9+ and 6_ ()3.14|) one can rewrite the quasilocal mass ()2.33|) 
as 

^ = f (1 + ^0+0-)- (3.17) 

so on any marginal two-sphere the quasilocal mass equals R/2. 

Now, which kind of global structure do Eqs. ()3.3|1 - ()3.4|1 allow for? We start by 
rewriting these equations in the integral forms 

p 

^'iP) = [ sm(20)f/p, (3.18) 





p 



2 A 

Q'ip) = I R'dp, (3.19) 





p 



R\p) = l-V j R^'^dp. (3.20) 



We see, that the first derivatives stay finite for any finite value of p as long as 
the metric functions R or Q don't vanish. If we assume the horizon Q{ph) = 

PH 

to be regular (and non-degenerate, Q'jj 7^ 0), we have to enforce / sin 20 dp = 



and therefore 

sin20H 

As Q is monotonically decreasing with p, pn is the only horizon along a m = const 
slice. 

Now the areal radius R of the two-spheres can behave in two essentially different 
ways: it can increase monotonically for all p or it can have an extremum at some 
Pe- 

In the first case, i?'(p) > Vp > 0, the solution exists for all finite p. From Eq. 
()3.5|) we see, that < i?' < 1 Vp > 0. Therefore R diverges like R = 0{p) for 
large p. From Eq. ()3.4|) follows, that Q' = 0{p) and therefore Q = O^p^). Eq. 
()3.3|1 shows, that goes to a constant at infinity. In this case the Killing horizon 
is a past inner trapping horizon (see Fig. 13.1k ) 

In the second case R'{pe) = 0, as R"{p) < for all p and -R'(p) < for p > pE, R 
goes to zero at some finite ps > pu (remember that we excluded the possibility 
of a second zero of R in the static region in Sec. I3.1.2|) . Again time evolution 
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beyond the horizon exists for all p < ps- However in the limit p ^ ps spacetime 
singularity occurs. This can be seen e.g. by investigating the behavior of the 
Kretschmann invariant 

+ 8Qi?' {R'Q'R" - R!) + 4Q2 {R!^ + 2R^B!'^) ) . (3.22) 

Since by assumption Q(p), Q\p) and R\p) are negative in a neighborhood of ps, 
all terms are non-negative and clearly not all denominators vanish in the limit 
p ^ Ps- Therefore the Kretschmann invariant blows up like 1/R^ in the limit 
P Ps- 

Concerning the character of the Killing horizon, two possibilities arise depending 
on whether R attains its maximum in the dynamic or in the static region. If 
Pe > Ph then the Killing horizon again is a past inner trapping horizon (See 
Fig. 13.11 b). If on the other hand pe < Ph, then the inward and outward 
direction interchange in the static region and the horizon therefore is a future 
outer trapping horizon and corresponds thus to a black hole horizon (See Fig. 

mt). 

The three possibilities are summarized in Figs. 13.1b - c. Fig. 13.1k shows the 
asymptotically de Sitter spacetimes. The areal radius R is increasing along a 
u = const slice from zero at the origin to infinity at X"*". A cosmological horizon 
separates the static region from an expanding dynamic region. 

Fig. 13. li b shows the situation, where R develops an extremum in the dynamic re- 
gion. Beyond the cosmological horizon spacetime initially expands until it reaches 
its maximal spatial extension at time pe and then recoUapses to a singularity at 

PS 

finite proper time ts = J dp/Q{p). 

PH 

Fig. 13.11 c finally shows a spacetime, where i?' = in the static region. At 
this hypersurface the inward and outward directions interchange. Therefore the 
collapsing dynamic region is enclosed by the static region. The separating horizon 
can be interpreted as a black hole horizon. 

We close this section by giving an upper bound for R{pe) for re-collapsing uni- 
verses (Fig. 13.1b ). As already mentioned the dynamic region of our solutions 
corresponds to Kantowski-Sachs universes. Therefore we can follow the work 
of Moniz [SHj on the cosmic no-hair conjecture of Kantowski-Sachs models. We 
re-investigate Eq. ()3.6|) . Remember that this equation is l/R"^ times the (p com- 
ponent of the Einstein equations. As the coordinate p plays the role of a time 
coordinate in the dynamic region, this equation is the Hamiltonian constraint for 
the time evolution problem. We rewrite Eq. ()3.fij) (devided by R^) as 
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Figure 3.1: Penrose diagrams for the three different global structures of solutions. 
Any point in the diagram (except where R — Q) corresponds to a two-sphere. The 
description is given in the text. 
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In terms of geometrical quantities on the spacelike hyp ersurf aces p = const, the 
terms on the left hand side are ^{K"^ — K''^ Kij) , Kij being the extrinsic curvature 
and K the trace of Kij. The last term on the right hand side is the energy 
density T^^n^n" and = \^^^'R' with ^'^'>1Z being the three scalar curvature of 
the slices. 

Note that the left hand side can only change sign if R' changes sign. So if at some 
instant of time po the universe is expanding, R'{po) > 0, then the left hand side 
is positive. Furthermore, as the energy densitiy is positive, we can give a lower 
bound for the left hand side, namely 

Suppose now that in addition at po the scalar curvature is smaller than 2 A, 
i.e. -R(po) > then the right hand side is positive. As -R'(po) is positive 

initially the lower bound on the right hand side increases away from zero, therefore 
making a change in sign of R' impossible. We can conclude from this, that if 
R'ipo) > and R{po) > l/y/X at some time po, the universe will expand for ever 
and approach de Sitter space for late times. On the other hand, a recoUapsing 
universe must have R{p) < 1/a/A for all times pn < p < Ps- 

3.2 Numerical Construction of Static Solutions 

The problem of constructing static solutions is given by the boundary value prob- 
lem Eqs. fj3.3|) - ()3.5p together with the regularity conditions at the origin ()2.28|) 
and at the horizon ()2.4(ij) . In order to determine the global structure of such a 
solution the data at the horizon (defined by the solution to the above boundary 
value problem) are used as initial data for the time evolution problem in the 
dynamic region. 

We solve the boundary value problem numerically using a standard shooting 
and matching method (see Appendix |XJ provided by routine d02agf of the NAG 
library jHOj- The numerical integration in the dynamic region is performed by 
routine d02cbf of the NAG library. 

In all our numerical calculations we set A = 3. 

3.2.1 Phenomenology of Static Solutions for A > 

In order to investigate the dependence of the spectrum of solutions on the coupling 
constant ?7, we start at ?7 = and follow a solution to larger values of 77. 
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For ?7 = 0, the vacuum Einstein equations with positive cosmo logical constant 
yield de Sitter spacetime, 



R{p) = p, Q{p) = 1 - 




(3.25) 



It was already shown in that the remaining boundary value problem for 
the field on fixed background admits a discrete one-parameter family of regular 
solutions. In the static region these solution oscillate around 7r/2. The number 
of these oscillations n can be used to parameterize the family. Energy increases 
with the oscillation number and converges to the energy of the "singular solution" 
= 7r/2 (see Sec. I3.1.3P from below. Outside the cosmological horizon at p = 
^/3/A the field goes to a constant at infinity (See Fig. 13. 2|) . 



3 




log(p/P//) 



Figure 3.2: The first three excitations 0„ on fixed de Sitter background. Within 
the static region the field oscillates n times around n/2, outside the horizon, it 
goes to a constant. 

We concentrated on the first three excitations and investigated their behavior as 
the coupling r] was turned on. We find that each member of these excitations 
exists up to a maximal coupling constant rjmaxin). The limit rj — > rjmax needs 
some care and will be discussed in Sec. 13.2.21 In the following we describe the 
properties of the solutions for < < rjmax- 

• For small ?7, < ?7 < rjcritiji) spacetime is asymptotically de Sitter as 
sketched in Fig. 13.1k and described in the accompanying text. The field 
behaves similar as in the uncoupled case. 
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• For Tj = r]crit{n) the areal radius R does not diverge like p but rather goes 
to a constant at infinity. For even bigger couphngs, rjcruin) < r] < rj^{n) 
the maximum of R occurs at earher and earher times pE (see Fig. 13. 3|) and 
time evolution ends in a singularity, as described in Fig. 13.1b . 

• Air] = 7]^ {n) the hypersurface, where R is extremal merges with the horizon. 
For bigger r], rj^{n) < r] < rjmaxin), R attains its maximum in the static 
region. The situation is as in Fig. 13.1b . 

Figure ()3.4|1 shows the areal radius R{p) for the first excitation for t] close to rjcrit- 
We have shown in Sec. l3.1.41 that re-collapsing universes must have R{p) < l/y/A. 
for all pH < p < Ps- Fig- 13.41 shows, that R{pe) for the re-collapsing universe 
comes close to the upper limit. This suggests that the upper bound for R is 
attained in the limit i] — > r/crit at infinity: R{oo; rjcrit) = 



R{p) 1 - 




1e-05 1e-04 0.001 0.01 0.1 1 10 100 1000 10000 100000 

p 



Figure 3.3: The metric function R{p) for the first excitation at couplings 
0.470366 < Tj < 0.470373. At rjcrit the global structure of spacetime changes 
from asymptotically de Sitter to a spacetime, that ends in a singularity. The 
vertical line denotes the cosmo logical horizon at p = 0.88761. 

The critical values of the coupling constants, ricrit{n),ri^:{n) and rjmaxin) depend 
on the excitation number ra, as can be seen from table IHTTl 



3.2.2 The Limit rj rj^ax 

Recall from Sec. 13.1.21 that the cosmological constant A sets the length scale in 
Eqs. fl3.3|l - ()3.5j) and that it can be eliminated from these equations, by intro- 
ducing the dimensionless quantities p = ^/Ap and R = ^/AR. This corresponds 
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Figure 3.4: The areal radius R{p) for the first excitation in the dynamic region. 
Plotted are the solutions at = 0.4703702 and t] = 0.470370903, close to r/c«t(l). 
The dashed horizontal line denotes 1 / -\/A (A was set to 3 for the numerical calcu- 
lation). As can be seen, the areal radius of the re-collapsing universe comes close 
to the upper bound 1/VA. We argue that therefore and due to the continuous 
dependence of solutions on the coupling rj, the upper bound is attained at infinity 
for 7] = T]crit- R{oo; r]crit) = 1/VA- 

to measuring all quantities that have dimension of length, as e.g. the energy E, 
the coordinate distance of the horizon pn from the origin, the radial geometrical 
distance of the horizon dn from the origin, the areal radius Rh of the horizon, 
and 1/0'(O), in units of I/a/A. We find that all parameters, that have dimen- 
sion of length go to zero in the limit t] f]max when measured with respect to 
this length scale. This indicates that is not the appropriate length scale 

for taking this limit. We therefore switch to the alternative viewpoint with pn 
as our length scale, and we fix pn = 1. In this setup A depends on rj and the 
excitation index n and goes to zero in the limit 1] — > rjmax- The parameters E, dn 
and 1/0'(O) attain finite values when measured in units of ph-, whereas Rh/ Ph 
goes to zero. (See Fig. 13. 5|) This strongly suggests, that there exists a solution 
with rj = rjmax which obeys Eqs. fl3.3|) - ()3.5|) with A = and has two centers of 
symmetry. In particular this means that the static region of this solution has no 
boundary, since any t = const slice has topology S'^. 

Furthermore, as can be seen from Fig. 13. (j^ the dimensionless parameter (plpn) 
for the first excitation tends to vr, and R'{ph) tends to —1 in the limit rj — > rjmax- 
As will be shown below, A = implies Q = 1. The limiting solution with A = 
will therefore satisfy the regularity conditions ()2.28|) and ()2.46|) not only at the 
axis p = but also at the second zero of R, which means that such a solution 
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n 




'/* 




1 


0.47037 


0.533 


1.0 


2 


0.41981 


0.48 


0.74255 


3 


0.41932 


0.474 


0.64931 


4 


0.42606 


0.4765 


0.60260 



Table 3.1: The critical values of the coupling constant for the first four excitations. 
The values of 77^, (n) were read off by eye and therefore are not as accurate as the 
other values. The values of rimax{n) are determined as described in Sec. I3.2.2| in 
particular the value of rj^axi^) = 1 is exact. While rimaxin) decreases with the 
excitation number n the other critical values do not seem to share this behaviour. 

is globally regular with two (regular) centers of spherical symmetry. In fact, for 
the first excitation this limiting solution is just the static Einstein universe ()3.9|1 . 
which was already given in closed form in Sec. 13.1.31 

These observations allow one to determine the maximal value of the coupling 
constant rjmaxin) not as a limiting procedure rj — >• rimax, but rather by solving 
the boundary value problem Eqs. ()3.3j) - ()3.5j) with A = and with boundary 
conditions, that correspond to two regular centers of symmetry. 

For A = Eq. ()3.4|) can be solved immediately to give R^Q' = const. According 
to the regularity conditions at the axis ()2.28|1 the constant has to vanish, which 
means that Q' = and therefore Q = 1. The remaining system of equations is: 

= sin(20) (3.26) 

R" = -r]R(j)'^, (3.27) 

and 

2ri sin^ - riR^^'^ + R'^ - 1 = 0. (3.28) 

Note, that this system of ODEs is scale invariant, that is any solution R{p), 4>{p) 
leads via rescaling to the one parameter family of solutions given by aR{ap), (p{ap). 
Keeping this in mind, we can fix the scale arbitrarily, e.g. in setting the first 
derivative of the field equal to one at the origin: 0'(p = 0) = 1. Thereby 
any solution, that is regular at the origin, is determined entirely by the value of 
the couphng constant rj. Regularity conditions at the second "pole" R{pp) = 
are the same as at the origin, except that either tends to vr, if its excitation 
number is odd, or to 0, if it has even excitation number. This can be inferred 
from 7r/2 < (f)H < vr for n odd and < (pn < 7r/2 for n even for all rj < rjmax, 
which is true for 17 = and according to ()3.21|) there cannot be any 77 < rj^ax 
where a solution with n > 1 has 0jy = O, 7ror7r/2 (apart from the "singular solu- 
tion" = 7r/2)^. Note that this corresponds to all odd solutions having winding 
number 1, whereas even solutions are in the topologically trivial sector. 

''if (pH = 0, TT or tt/2 then = (f)'^ = and therefore all higher derivatives vanish 
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These regularity conditions together with the invariance of the equations under 
reflection at the location of the maximal two-sphere R\pe) = 0, causes globally 
regular solutions R{p) to be symmetric around pe whereas (j){p) — 7r/2 is either 
antisymmetric for n odd or symmetric for n even. For (j) symmetric the formal 
power series expansions of -R(p) and 0(p) around p = Pe gives 



Rip) 



R{pj,) + 0{{p ~ peY 



arcsm 



2v/l- 1/277 {p-pEf 



R{Pe) 

and for </) — 7r/2 antisymmetric we get 



2! 



+ 0((p-p^r), (3.29) 



R{p) 



v(i^'{pEy 

IT 



2! 



:^ + ^'{pj,){p-pj,) + 0{{p-pEf). 



(3.30) 



In order to solve the system ()3.26|) . ()3.27|) we again use the shooting and matching 
method on the interval [origin, pe] using the above Taylor series expansions to 
determine the boundary conditions at p = pe- Shooting parameters are now 
Pe,(P'{pe) and r] for odd solutions and Pe,R{pe) and i] for even solutions. The 
results are displayed in Table |S21 



n 




Pp = dp 




E/2nf^dp 


1 


1 


TT 


37r/2 


3/2 


2 


0.74255 


6.74225 


11.78039 


1.74724 


3 


0.64931 


12.10140 


22.43662 


1.85405 


4 


0.60260 


19.63717 


37.47302 


1.90827 



Table 3.2: Results for the first three excitations for A = 0. Since Q = 1 the 
coordinate distance pp of the two regular "poles" equals the radial geometrical 
distance dp. The energy density pp and energy E are given in units where 
0'(O) = 1. The ratio E/dp can be compared to the results for solutions with 
A > and represents the limit r] — > rimax for those solutions . 

It is clear from ()3.29j) and ()3.30|) . that regular solutions for A = can only exist 
if ?7 > 1/2. Assuming now that our numerical observations concerning the first 
few excitations extend to higher excitations, we argue as follows: since every 
"branch" of the "A > solutions" persists up to a maximal value of t], which 
can be computed by solving the boundary value problem ()3.26|1 . (I3.27|1 together 
with regularity conditions at the two "poles" - which implies rj > 1/2 - and since 
we know, that in the limit r/ — > there exists an infinite number of excitations 
we conclude that this whole family of solutions with A > persists up to 
some maximal value fjmaxi'^), which is greater than 1/2. In other words, for any 
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Figure 3.5: Some parameters of the first excitation - measured in units of ph - 
as functions of the couphng constant. Except for Rh the parameters attain finite 
values in the hmit rj — > rjmax, when measured with respect to this unit. 




V 



Figure 3.6: The parameters (I){ph) and R'{ph) for the first excitation as functions 
of T). In the hmit rj —>■ rj^ax <P{ph) tends to tt and R'{ph) goes to —1, which are 
necessary conditions for a second regular center of spherical symmetry at pH- 



43 



7] < 1/2 there exists a countably infinite family of solutions with A > 0, whereas 
for T] > 1/2 our numerical analysis shows, that only a finite number of solutions 
exists. (See Table . 



3.2.3 Stability 

In order to analyze the stability properties of the above described static solutions, 
small perturbations of the metric functions and the field are considered. We set 

Q{t,p) = Qn{p) + SQ{t,p), 
R{t,p) = R^{p) + 6R{t,p), 

= 0„(p) + (50(t,p). (3.31) 

Qnip), Rn{p) and (pnip) denote the n-th static excitation. The perturbations 
SQ{t,p), 6R(t,p) and 6(f){t,p) are considered to be small such that the equations 
can be linearized in these quantities. 

Inserting the ansatz 1)3.311) into Eqs. ()2.51|1 - ()2.54|1 . making use of the static 
equations and linearizing in the perturbations gives a coupled system of linear 
PDEs for the perturbations. The momentum constraint ()2.54|) contains only first 
time derivatives, which enter each term linearly, so this equation can be integrated 
with respect to time to give 

QM6R{t, p) 2Q4t, p)SR'{t, p) 2r]Q4p)R4p)<l)M6<j){t, p) 



6Q{t,p) 



Kip) Kip) R'nip) 

(3.32) 



The constant of integration, a function g{p) is determined by the Hamiltonian 
constraint (making use of the other equations) to be 

const 

Rn{p) R'nip) 

so for perturbations SQ{t, p) regular at the origin, the constant has to vanish and 
therefore g{p) =0. 

Eqs. (ITH^ and (IT^ then give 

6Rit, p) = -Qlip) {6R"it, p) + v€ip)mt, P) + 2(5Rn{p)(t>M5ct>'{t, p)) 

(3.34) 

and 

-^^^*^ + (i?^(p)Q^(p)^0'(t,p))' + 

+ {^'^§J^ iQM6R{t,p) - 2Qnip)SR'{t,p) - 2r^Qnip)Rnipmpmt,P)) 
+ {2Rnip)QnipWnip)SRit,p))' = 2cos(20„(p))50(t,p). (3 
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Now it is important to note, that the form ()2.17j) still contains a certain gauge 
freedom. Therefore the perturbations ()3.31|) may - in addition to physical modes 
- also represent pure gauge modes. We investigate this, by considering small 
coordinate transformations 

t t = t + ex\t,p) (3.36) 
P ^ p = p + ex'it,p). (3.37) 

As described in Sec. 14.3.21 the perturbations ()3.3H) then transform up to order e 
according to the Lie derivative along x of the "background solution". In detail 
we have 

5Q = SQ-eiQ'^x' + ^QnX't) 
6R = SR-eR'^x" 

H = 50 (3.38) 
where x is subject to the following conditions: 

^X't = -QnXl (3.39) 
X't = -X'p. (3.40) 

Condition ()3.39|) comes from the fact, that the shift is zero in this gauge and 
therefore C^{gn)pt = 0, and condition ()3.40|1 arises from the fact, that 6gtt = 
QlSgpp and therefore C^{gn)tt = QlC^{gn)pp- Combining Eqs. (!3.39jl and (!3.4(ljl 
yields a wave equation for x^: 

Qlx'pp - X% = 0. (3.41) 

Gauge transformations that respect the choice ()3.3H) therefore are determined by 
the single function x'', which is subject to Eq. 1)3.411) . 

Pure gauge modes are perturbations, that can be removed by a coordinate trans- 
formation, i.e. they have to be of the form 

5Q - e{Q',x' + 2QnX't) 
6R = eKx' 

50 = e^W, (3.42) 

with x^ and x* subject to ()3.39I3.40|) . Clearly pure gauge modes satisfy the 
perturbation equations ()3.34|) - ()3.35p . 

On the other hand, one can try to combine the perturbations in such a way, that 
the combination is invariant under a coordinate transformation. One can show, 
that such a quantity has to be of the form 

e(p,t) = aip)iR',ip) 6fip,t) - ap) 5R{p,t)), (3.43) 
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where a{p) is an arbitrary function of p, which we choose to be l/i?„(p). Using 
the definition for ^ and Eqs. ()3.34|) . ()3.35|) we obtain a single pulsation equation 
for the gauge-invariant quantity ^: 

- ^"^^^ P) + Qn{p) Rn{p? R'Ap) eit, p) + (2 QM Rn{pf ^ip)' + 
Qn{p) V 

+ Rnip)' QM Kip) + 4 Quip) Rn{p) K{pf) i\t. P) + 

r, Rn{pf QM - 2 cos(2 0„(p)) <(p) + 

+ 2 g„(p) i?„(p)2 0;(p)2 K(p) + 2 gn(p) - 

- i?„(p) (2r7sin(20„(p))0'„(p)-Q;(p)i?;(p)2))e(t,p) = O. (3.44) 

As the coefficients do not depend on time, we can work with Fourier modes 

e(t,p) = e*'^*l/(p), (3.45) 
which turns Eq. ()3.44|1 into a linear second order ODE 

Q„(p)2 R^{pf RM y"ip) + (2 V Quip)' Rn{pf Up? + 

+ Qn{p) Rn{pf Q'nip) K{p) + 4 Qn{pf Rn{p) Kipf) V {p) + 

+ {-2r^Qn{p)Rn{p) Sm{2<Pn{p))<P'n{p)+VQn{p)Rn{pf<P'n{pfQ'n{p)- 

- 2 cos(2 0„(p)) g„(p) K(P) + 2 g„(p)2 Rn{pf 4^M^ K{p) + 

+ Qn{p)Rn{p)Q'n{p)K{pf'r1Q^{pfK{pf)y{p) = ~CT^ Rn{pf K{p)y{p)- 

(3.46) 

Again this equation has regular singular points at the origin where -R(O) = and 
at the horizon Q{ph) = 0. The corresponding regularity conditions for y{p) are 



l/'(0) = 0, yip)r^ipH-pr, with a = -^^, (3.47) 
for negative cr^. 

If the background solution has i?^ > for all p < pn, then Eq. 13.461 together 
with the boundary conditions constitutes an eigenvalue problem with eigenvalue 
and eigen vector y{p). As we are interested in unstable modes, we look for 
negative eigenvalues a^. 

If on the other hand the background geometry contains a maximal two-sphere 
R'{Pe) = within the static region pe < Ph, then Eq. ()3.46|) has an addi- 
tional singular point within the range of integration. The behavior of the two 
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independent solutions in the vicinity of this singular point is 



oo 



1/1 (p) 



{p- Psf'^anip- Pe) 



n 



n=0 



oo 



1/2 (P) 



(3.48) 



so the general solution stays regular near p = pe- Nevertheless, the coefficients 
of the ffist and zeroth derivative of y{p) in Eq. ()3.46|1 are singular, which renders 
the standard numerical shooting and matching methods impossible. 

Nevertheless we tried to solve this problem, using a standard relaxation method 
(routine d02raf of the NAG library [BOJ) on one hand. On the other hand, we 
discretized Eq. fl3.46|) by hand, thereby turning the Sturm Liouville eigen value 
problem into an algebraic eigen value problem. We only present preliminary 
results here: for 77 = stability was already analyzed in 0. It turned out, that 
the n-th static excitation has n unstable modes. Our preliminary investigation 
of the first three excitations in the coupled case gives, that these solutions keep 
their unstable modes when gravity is turned on until in the limit of maximal 
coupling one negative eigenvalue crosses zero. The occurrence of one "indifferent 
mode" in the hmit 77 = rimax is due to the fact, that there the equations are 
scale invariant, i.e. for each coupling rjmaxin) we have a one parameter family of 
solutions (-Ra(p) = Ai?„(Ap), (f>\{p) = 0„(A)(p)). For A = 1 + e we can write 



Up) 



Rn{p)+e{Rn{p)+R'o{p)p) 
<Pn{p) + e0'„(p)p. 



(3.49) 
(3.50) 



which just corresponds to an "indifferent" mode with cr = 0. 
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Chapter 4 



Self-similar Solutions of the 
Self-gravitating a Model 

As self-similar solutions usually govern type II critical phenomena (see Chap.Ej), 
it is of great advantage to know about their properties like existence and stability 
from a "direct construction". The "direct construction" - in contrast to a con- 
struction that uses time evolutions of near critical data (see Chapter - profits 
from the symmetry, which is imposed on the equations. This way the problem 
of constructing self-similar solutions reduces to a boundary value problem for 
ODEs in the case of continuously self-similar solutions (see Sec. I4.2|l and to a 
time-periodic boundary value problem for PDEs in case of discrete self-similar 
solutions (see Sec. 14. 4p . CSS solutions of the a model have already been con- 
structed by Bizon and Bizon and Wasserman [TTj . We reproduce their results 
here. Stability is analyzed via the usual method for CSS solutions fSec. l4.3.T|) . In 
addition J. Thornburg [Slj proposed and carried out a method, that is based on 
discretization of the field and uses the full (nonlinear) field equations (Sec. l4.3.T?jl . 

We summarize our results in Table I^^Tl In Sec. 14. 5.11 finally we compare the CSS 
(first excitation) and DSS solution in a certain range of couplings, where both 
solutions exist. Our observations lead us to conjecture, that the DSS solution 
bifurcates from the first CSS excitation in a homoclinic loop bifurcation at ?7 ~ 
0.17. To our knowledge up to now such a bifurcation has not been observed in 
the context of self-similar solutions to self-gravitating matter fields. 

4.1 Continuous and Discrete Self-Similarity 

A spacetime {M,g) is said to be discretely self-similar (DSS) if it admits a dif- 
feomorphism $a '■ M M, which leaves the metric invariant up to a constant 
scale factor, that is 

i^l9)\p = e'''g\, ypEM, (4.1) 
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where A is a real constant. 

A spacetime, that admits a one-parameter family of such diffeomorphisms, parametrized 
by A and with $o = idu, is called continuously self similar (CSS). The gener- 
ating vector field ^ = ^ '^aIa=o ^ special case of a conformal Killing vector 
field. It obeys the conformal Killing equation 

9 



with a constant factor in front of the metric at the right hand side. A vector field 
satisfying ()4.2j) is called homothetic. 

At any point p G M we can compare the original metric g\p and the metric g that 
results from pulling back 5'l$^(p) to p, g\p = ($A5')lp = g\p. As g and g only 
differ by a constant rescaling both metrics yield the same covariant derivative, 
Vg = Vg. Therefore the Riemann as well as the Ricci tensors are identical, 
'^%ui9)\p = T^''^rui9)\p and TZf,^{g)\p = n^^{g)\^. The Ricci scalar scales as the 
inverse metric, 

n9)l = m,M\p = e-^V^7^,.(^7)|^ = e-2^7^(^7)|^, (4.3) 
and therefore the Einstein tensors for both metrics agree: 

G,M\p=G,Ag\. (4.4) 



Now the Riemann, Ricci and Einstein tensor and the scalar curvature (summa- 
rized as T) behave under the pull-back of a general diffeomorphism f : M ^ M 
as 



ifTmp 


= nrg)\p. 


(4.5) 


For a DSS spacetime we therefore have 














= TZ 1 














(4.6) 



For a CSS spacetime the above considerations directly yield 



L(]z = {c^g-'rn^, = -2n, 



= I-oa(^i^-^^*-)*^)I.) = 

= l^ox(^l.-^""^l.)=2 4- (4-2) 



(4.7) 
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4.1.1 Adapted Coordinates 



In order to simplify the discussion, we introduce coordinates, that are adapted 
to the symmetry (see e.g. 

In order to formally construct such a coordinate system for a DSS spacetime, we 
choose a hypersurface S, with S fl $a(S) = 0, provide it with coordinates (z*), 
i = 1, 2, 3 and label it with tq. (Up to now this hypersurface can be spacelike, null 
or timelike). We apply the diffeomorphism to S, label $a(S) with tq — A and 
choose the coordinates in this hypersurface such that $a(to,-2*) = (tq — A,^*)- 
Next we choose coordinates (r, z'') in between these two hypersurfaces, with tq — 
A < T < Tq and their restriction to S and $a(S) being (ro,z*), respectively 
(to — A,z^). Then we use the diffeomorphism to copy this coordinate patch to 
the other regions of spacetime. Of course this construction is very far from being 
unique. 

Per construction the diffeomorphism maps a point p = (r, z^) to the point 
$a(p) = (t — a, z^). So the Jacobian is the identity, = 6^. 

Using the definition of the pull-back as well as Eq. ()4.1|) in these coordinates we 
obtain for p = (r, 2;*) 

^ e'%Ar,z'). (4.8) 

As this is valid for any point p & M we can conclude that the metric is conformal 
to a metric g, which is periodic in the coordinate r, with the conformal factor 
being an exponential in r, 

gt^A^, z') = e~^^ g^^T, z') with g^^r + A, z') = g^^T, z'). (4.9) 

For a CSS spacetime again we choose a hypersurface S, with En$A(S) = VA, so 
the homothetic Killing vector field is transversal to S. We choose coordinates (z*) 
on this hypersurface and transport them across spacetime via the one-parameter 
family of diffeomorphisms. We parameterize the orbits 7^ of the homothetic KVF 
^ with — r, so ^ = —dr- The freedom in the construction in this case consists of 
choosing the hypersurface and applying diffeomorphisms within S. 

In analogy to ()4.8|) we have 

{^i9)f,u = -drgf,u '29f,u, (4.10) 

and therefore 

g,AT,z')=e-'^g,Az'). (4.11) 
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4.1.2 Self- Similarity in Spherical Symmetry 



The fact that the above introduced adapted coordinates are not unique is no 
draw back for our purposes. What we are interested in is to find out whether our 
numerically evolved spacetimes contain self-similar regions. In order to do so, 
we need to know, what self-similarity looks like in the coordinates our code uses, 
namely the Bondi-like coordinates in spherical symmetry, defined in Sec. fl2.2j) ). 

ds^ = -e^'^^'^^^Uui-iu, r)du + 2dr) + r^dfil uu) 
r 

Assuming that ()4.12|) describes a self-similar spacetime, we seek for a coordinate 
transformation (u, r) — >■ (t(m, r), r)), such that the resulting metric is of the 
form e~'^'^g^u, where g^u is periodic in r for a DSS spacetime and independent of 
T in the case of a CSS spacetime. In the following a "~" means, that the function 
is periodic in r with period A. 

The first observation is, that the coordinate transformation does not involve the 
angles 9 and ip. So the "S^"-part of the metric is unchanged. We immediately 
get 

r{T,z) = e-^R{T,z) for DSS, (4.13) 
r{T,z) = e-^R{z) for CSS. (4.14) 

This means that the diffeomorphism maps an r = const hypersurface to the 
hypersurface e^r = const. 

In the following we exploit the relations ()4.1|) and ()4.2j) with respect to the Bondi 
coordinates, in order to find the behavior of m = const hypersurfaces under the 
diffeomorphism as well as determining the r dependence of the metric functions 
P and ^. 

For the DSS spacetime we write = ($^(-u, r), $^(m, r)) = ($^(m, r), e^r), 

the last equality following from ()4.13|1 . Therefore d(^\/du = 0, d^^^/dr = e^. 
We first examine the rr component of ()4.H) : 

{^lg)rr = e^^grr = (4.15) 

and so 

-gf [-dfSuu + 2e^gur] = 0. (4.16) 

This formula states, that the null vector V^m, which is tangent to the outgoing 
null geodesies generating the u = const hypersurfaces, is mapped again to a 
null vector. If the first factor vanishes, = 0, the push forward of 'V'^u\p 
is parallel to V'^m|$^(p). If the expression in the parentheses vanishes the push 
forward would be parallel to the ingoing null geodesic vector V^f |$^(p). Here we 
are interested in those diffeomorphisms that are connected to the identity map. 



51 



Therefore we want the first factor to vanish, so = and $^ = This 
shows, that u = const hypersurfaces are mapped to m = const hypersurfaces. 

We invoke now e.g. the {u, r) component of ()4.1|) in order to get more information 
on 



-^y^'guri^^ip)) = e^^gUp)- (4.17) 

Consider now a point p at the origin, i.e. p = (m, r = 0). As the origin is mapped 
to itself (which follows from ()4.13|l ) the image of the point p is $a(p) = {^\{u), 0). 
As described in Sec. (12.2.11) the components of the metric with respect to Bondi- 
like coordinates are fixed at the origin, due to regularity at the center as well as 
the choice of retarded time being proper time at the origin. In particular we have 
P{u,r = 0) = 0. Inserting this into Eq. ()4.17j] we obtain 

Integrating gives 

(^l(u) = e^u + const. (4.19) 
Reinserting this into ()4.17|1 we get 

emr-A,z) ^^2(3ir,z)^ (4.20) 

and for the {u, u) component of ()4.H1 we have 

i^*Ag)uuL = e^^gn 



-^j guui^Aip)) = e'^'gUp) 
e2/3(-A,.);^(^_A,^) = e2^(^'^)^(r,z). (4.21) 

From ()4.2()|1 and ()4.21|1 we now see, that both metric functions are periodic func- 
tions in r 

(5{u,r) = l3{T{u,r),z{u,r)), 
V V 

— (M,r) = — (r(M,r),z(M,r)). (4.22) 

Note, that this relation is valid for any set of adapted coordinates constructed as 
described in Sec. fj4.1.1|) . 
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For a CSS spacetime we use analogous arguments. We start by writing the 
homothetic Killing vector as 

-drl=i = r)du + r)dr. (4.23) 

Of course we have = —drU and = —drV (both partial derivatives taken at 
constant 2;*). ()4.14p immediately gives the r component of 1^, 

C = -drr{T,z)\^ = r{T,z). (4.24) 
As above we examine the rr component of ()4.2|) . 

C^Qrr = e9rr,f. + 2 ^^r9^.r = 2 ^"^rQur = 0. (4.25) 

So we have = and ^" = ^"(m). Now the uu component of the homothetic 
Killing equation ()4.2|) reveals 

^i;9uu ^^9uu,fi ~l~ 2 ^^(jfj.u 2 gy^u 

e9uu,f. + 2 Qguu = 2 gun. (4.26) 

Again we use the fact that the origin is an orbit of the homothetic Killing vector, 
so ^1^=0 = du\r=o, and guu\r=o = 1- So (I4.26|) gives 

= 1 ^ ^"(u) =u + const. (4.27) 

Therefore Eq. ()4.26|) gives 

dr9uu = 0, (4.28) 
and the (ur) component of ()4.2|) yields 

^(,9ur ^^9ur,fi ~l~ ^^9 fir ~l~ ^^9ufj, 2 

^^9ur,fi ~l~ '^9ur 2 5'n,r 

= 0. (4.29) 

So in analogy to ()4.22|1 we have for a CSS spacetime 

f3{xn = 

V V 

-^{x^) = -^{z). (4.30) 

We now explicitly write down a set of adapted coordinates for DSS and CSS 
spacetimes, which will be used later. 

As we have seen, the diffeomorphisms map a. u = const hypersurface to another 
u = const, different form the first. Therefore the u = const hypersurfaces are 
valid candidates for r = const hypersurfaces, i.e. we may set r = t{u). 
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For a DSS spacetime the function u{t) has to obey 

e\{T) + const = u{t - A). (4.31) 
The simplest function u{r) satisfying ()4.3H1 is given by 

u* -u = e-^, (4.32) 

where we have set u* = const e^. 
For the CSS spacetime we integrate 

du 

- — = C = u + const, (4.33) 

(Jj I 

which gives 

\n{\u + const\) = -T + a, (4.34) 

Ci - being a simple shift in r - can be chosen arbitrarily, and so we set ci = 0. If 
we again introduce u* = —const, we have 

u* -u = e-\ (4.35) 

the left hand side coming from resolving the absolute value in ()4.34j) for u < u*, 
which is the region we will be interested in in critical collapse situations. 

Finally we parameterize the r = const hypersurfaces as follows: for a DSS space- 
time we set 

r{T,z) = e-^R{T,z) = e-^ z C{r) with C(^ + A) = C(r), (4.36) 

where we require ({r) > ({r) for all r, such that drvir, z) < for all r. There 
is still a coordinate freedom contained in C(''")- We keep this in order to describe 
the null hypersurface, which is mapped to itself via the diffeomorphism (the so 
called self-similarity horizon) with z = const. The resulting condition on ({t) 
will be given below. 

For a CSS spacetime we set 

r = e-^R{z) = e-^z. (4.37) 

Summarizing, we chose the following adapted coordinates for a DSS spacetime 
t{u) = — hi{u* — u) u{t) = u* — e^^ 

ziu, r) = r{r, z) = e-^zdr), (4.38) 



54 



and for a CSS spacetime 

t{u) = — ln(M* — u) u{t) = u* — e""^ 



{u* — u) 



(4.39) 



By construction, the hypersurfaces r = const are null. For m — — oo we have 
r — oo, whereas for u u*, r — +00. 

The hypersurfaces z = const, along which the diffeomorphism acts all meet in the 
point {u = u*,r = 0). Clearly the adapted coordinates get singular there. Due to 
the symmetry, which squeezes the geometry into smaller and smaller spacetime 
regions, this point also has diverging curvature, as can easily be seen as follows: 
We examine the scalar curvature 71 in adapted coordinates. From ()4.6|) we have 



7^($A(p)) = e-2^7^(p) 
7^(r-A,2) = e-^^7^(r,^) (4.40) 



and therefore 

7^(r, z) = e^^n{T, z) with 11{t - A, 2) = 11{t, z). (4.41) 

Note that 'R. does in general not agree with the scalar curvature built from the 
metric introduced in Sec. ()4. 1.111 . although this would have a periodic r de- 
pendence as well. 

In particular, as 7?.(r, 2;) has a periodic r dependence, it is bounded for all r if it 
is bounded in one "segment" between E and $a(S). So moving along z = const 
we find the scalar curvature blowing up like e^"^ for t ^ 00, or u ^ u* (unless TZ 
vanishes identically of course). As the origin is a line of constant z, this blowup 
occurs at {u = u*,r = 0). We call the point {u*, 0) the culmination point, and u* 
the culmination time. 

Another point is worth to note. Consider the square of the vector field dr- From 
.38j) we have 

dr = {u* - u)d^ + r(-l + ^)dr (4.42) 



and so 



g{dr, dr) = {u* - ufguu + 2r(M* - m)(-1 + '^)guT = 

= _e-2-e2/^(^'^)(-(r, z) - 2<(r) + 2<(r))). (4.43) 
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At the origin (z = 0), dr is timelike, which is clear, as the diffeomorphism maps the 
origin onto itself. Furthermore, as mentioned after ()4.36|) we assume ({t) > ({t) 
for all r. 

Moving away from the origin z grows and goes to infinity for u ^ u* . As 
we assume ^ to be bounded, the expression in the parentheses in ()4.43|) will 
vanish at some value of z, which depends on r. We can now use the additional 
gauge freedom, which is contained in C,{t) in order to have dr getting null on a 
hypersurface of constant z, e.g. at z = 1. The resulting condition on (^(r) reads 

-(r,l)-2C(r) + 2C(r)') =0. (4.44) 
r / 

So this hypersurface, z = 1, is a null hypersurface, which is mapped onto itself 
via the diffeomorphism. It is called the past self- similarity horizon (SSH). In fact 
this self-similarity horizon is just the backwards light cone of the culmination 
point. 

There exists another null hypersurface which is mapped to itself by the diffeo- 
morphism. This is the future light cone of the culmination point, and is called 
the future SSH. As the culmination point corresponds to a spacetime singularity, 
the region beyond this future SSH is not determined by the solution given for 
u < u*. 

For a CSS spacetime the statements "periodic in r" have to be replaced by 
"independent of r". Furthermore with our choice of coordinates the location of 
the past SSH is given by 

-{zh)=2zh. (4.45) 
r 

Fig. 14.11 shows a schematic diagram of a self-similar spacetime. 

Note, that a self-similar spacetime is not asymptotically flat (unless spacetime as 
a whole is fiat). This can be seen by going to infinity on a spacelike hypersurface 
z = const outside the past SSH. As the metric functions are periodic in r (resp. 
constant) along these hypersurfaces, they do not fulfill the fall-off conditions 
required by asymptotically flatness. 

We close this section by writing down the line element for a DSS (CSS) spacetime 
with respect to the adapted coordinates ()4.38|) (respectively ()4.39|) ). 



_^2p(r,z) \^{t,z) - 2<(r) + 2<(r)| rfr^ - e^'^^^''hC{T)dTdz+ 
+ z^irfdn-"] , (4.46) 



and 



(^^css — e ^'^ 



-e^m J ^t,) _ 2z I dr' - e'^^^hdrdz + z'dn' 



r 



(4.47) 
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Figure 4.1: A schematic diagram of a self-similar spacetime. The adapted co- 
ordinates {t,z) defined in ()4.39|1 cover only the region u < u*. In this region 
the metric functions (3 and ^ are constant along z = const, thereby shrinking 
their profile to zero at the culmination point {u = u*,r = 0), where (unless in 
flat space) a spacetime singularity occurs. The region within the backwards light 
cone of the culmination point is the one of interest for critical collapse situations. 
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4.1.3 Implications for the matter field 

We turn now to the conditions that the symmetries ()4.1|) and ()4.2j) imply for the 
matter field. It is clear that the matter will have to share (at least part of) the 
symmetry imposed on the stress energy tensor via the Einstein equations, i.e. we 
have for a DSS spacetime 



($1T) 



T, 



and for a CSS spacetime 



C^T^u = 0. 



(4.48) 
(4.49) 



As defined in Sec. 12.11 the stress energy tensor of a harmonic map is 



Tfii/ (x 



(^V^X'^(x)V.X^(x) - ^g,,{x){g''^ix)V.X^{x)VrX''{x))^ G^ij(X(x)). 

(4.50) 

The pull-back of this tensor under $a is given by 

i^lT),^\^ = [V.X^($a(|>))V,X^($a(p)) - 

- ^^7„/5($a(p))(^^'($a(p))V,X-^($a(p))V5X^(<|.a(p))) 



Gab{X{<^a{p))). 
As the inverse metric transforms like 

(($A)*r'^^'l 



*a(p) 



4>a(p) ' 



(4.51) 



(4.52) 



we have 



^7^^($a(p))V,X^($a(p))V5X^(<I>a(p)) = 
e-^^-(p)^^V,X^($A(p))V,X^($A(p)). 



(4.53) 



The factor e cancels with e^^ coming from 



d<S>1 d<S> 

dxi^ dx 



^daisi^Aip)) in the second 



term of ()4.5H) . In order to simplify things, we contract ()4.5H1 with g^^{p) 



I 



dx 

and switch to adapted coordinates to get 



-'^r ,i V.X-^($A(p))V;3X^(«l>A(p))gAi.(X($A(p))), 



(4.54) 



(?'^^(r, z, a;) V^X^(r - A, z, a;) V,X«(r - A, z, u)Gab{X{t - A, z, u)) 



g^'^iT, z, cc.) V^X^(r, z, u;) V,X«(r, ^, uj)Gab{X{t, z, u)) 



(4.55) 
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where we abbreviated the angular coordinates 9 and Lp with u. So the self- 
gravitating harmonic map, leading to a discretely self-similar spacetime has to 
be either periodic or "anti-periodic" in the adapted time coordinate r, 

X^(a;^) = X^(r,2,cj) with X^(r - A, z, cj) = ±X^(r, z, tu), (4.56) 

and independent of r for a CSS spacetime, 

X\x^) = X\z,uj). (4.57) 

Concerning the sign in ()4.56p we adopt the following convention: if the metric 
functions are periodic with period A and the field is "anti-periodic" with respect 
to this period, then of course it is periodic with respect to twice the period, i.e. 
2A. In this case we say that the solution is DSS with period A = 2A, and has 
the additional symmetry /3(r -|- A/2,2;) = P{t,z) etc. and X^(r + A/2,z) = 
-X\r,z). 

Of course the conditions ()4.56|) and ()4.57|) are only necessary conditions for the 
existence of a regular self-similar spacetime^. 

Remember, that the self-gravitating nonlinear a model is a scale invariant theory. 
We stress that this scale invariance is a necessary condition for the existence of 
self-similar solutions. Consider for example vacuum with a cosmological constant 
A. 

G,.. + kg^u = 0. (4.58) 

Here the cosmological constant A introduces a length scale. Applying the pull- 
back $^ to Eq. we get using Eqs. (g^l) and dHH) 

($1(6- + Xg)^,\^ = G,,, + Ae2%,|p, (4.59) 

which does not satisfy ()4.58|) anymore. 

Another example of a model with a length scale would be the self-gravitating 
nonlinear a model with an additional potential in the Lagrangian, e.g. ViX^ = 
X'^X^GAsiX). For dimensional reasons this potential has to be multiplied by 
a constant of dimension {I/length^ , which breaks the scale invariance. As this 
potential term appears in the stress energy tensor multiplied by the metric g^^, 
it would get a factor e^^ under the action of the pull-back <I>^, which does not 
cancel. As this is the only term which transforms that way, again the Einstein 
equations cannot be satisfied. 

Nevertheless a model with a length scale can admit asymptotic self-similarity, i.e. 
it might display self-similarity at scales which are small compared to the length 
scale of the theory. This can be seen by writing the equations in adapted coordi- 
nates and neglecting the terms which contain a factor or any power thereof. 

^By regular in this context, we mean regular within the backwards light cone of the culmi- 
nation point 
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(As r — oo denotes the region close to the culmination point, spatial extensions 
are already very small). This concerns the terms that are tied to the length scale, 
therefore the remaining theory again is scale invariant. In this context such terms 
are called asymptotically irrelevant. An example is the Einstein- Yang-Mills sys- 
tem, which admits an asymptotically DSS solution at the (type II) threshold of 
black hole formation [21]. This solution also has been constructed directly (using 
the asymptotic symmetry) by Gundlach j34|| . 

4.2 Numerical Construction of CSS Solutions 

This section deals with the (numerical) construction of CSS solutions of the 
self-gravitating SU(2) cj-model in spherical symmetry using the hedgehog ansatz 
introduced in Sec. 12. 2. HI This problem has already be studied by Bizon |Hj for 
the simpler case of fixed Minkowski background {rj = 0), and by Bizon and 
Wasserman for the coupled case [TTj . 

In jH] a discrete one-parameter family of CSS solutions was constructed numer- 
ically, the existence was proven analytically and the stability properties of the 
solutions were given. In pTj this family of solutions was shown to persist up to a 
coupling constant rjmax = 0.5 by means of the numerical construction. Further- 
more, the analytic continuation beyond the past SSH was studied numerically, 
showing that for each member of the family there exists a critical value of the 
coupling ?7* beyond which the analytic continuation contains marginally trapped 
surfaces. 

In order to be able to compare the directly constructed CSS solution to critical 
solutions obtained by a bisection search of time evolved data, we re-did the (nu- 
merical) calculations of [Sj and fl], reproducing their results. Furthermore we 
studied the stability of the solutions for nonzero couplings, which has not been 
considered in [TTj . 

We will report here on both the numerical construction of the solutions and 
their analytic continuation beyond the past SSH (in order to make everything 
self-contained) as well as on their stability properties. 

4.2.1 The "CSS equations" 

We start by combining the symmetry requirement ()4.57j) with the hedgehog 
ansatz. For a CSS solution of the self-gravitating nonlinear a model we have 

with = (0, e, $) 

0(x'') = (f){z), e(x^) = e, ^x^) = ip. (4.60) 
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Transforming now Eqs. ()2.48j) . ()2.55j) and ()2.5(i|l to the adapted coordinates de- 
fined in ()4.39|) and dropping all derivatives with respect to r, we get 

(3' = ^z{<pr (4.61) 



and 



or 



r z \ r J 

{z^^ - 2^)0'^ + = sin(20)e2^ 



(-e^^ {sin(2 0) + ;z (-1 + 2r7sin(0)^) 0'} 



(4.62) 



(4.63) 



+ z {[--Az]<P'\], (4.64) 



where ' = d-^- 

The subsidiary Einstein equation "Euur" ()2.57|1 gives in addition 

(^)' = f 2^ - -V^')'- (4.65) 
r \ r J 

Eqs. fl4.62|) and ()4.65|) can be combined to give an algebraic relation for ^: 

V -62/3(1 -2r/sin2 0)+r^2z3(0')2 



-1 + rjz 



(4.66) 



As discussed already in Sec. 12. 2| regularity at the origin {z = 0) as well as the 
gauge choice for u requires 

0(0) = 0, P{0) = 0, -(0) = 1. (4.67) 

r 

The only free parameter, which determines the solution is 

0'(O) = b. (4.68) 

According to ()4.65|1 ^ is decreasing for 2; > and eventually equals 2z at some 
zh- This marks a singular point of the equations and corresponds physically to 
the past self-similarity horizon discussed in Sec. 14.1.21 

Eq. ()4.65|) shows, that at the horizon in addition to^ln = 2zh , we have \h = 
0. Eq. ()4.62|) then yields the following relation for jSn 

e^fe = ^'"^ (4.69) 

l-2r/sin2(0H)' ^ ^ 
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where [3h, 4>h denote the fields (3 and evaluated at the horizon zh- 

Furthermore following from Eq. ()4.fj4|l regularity at the horizon requires sin 2(f)H = 
0, which can be resolved to (pn = (mod vr) or (pn = 7r/2 (mod tt). The first case 
is impossible for the following reason: assume (pn = 0, so e^'^^ = 2zh which in 
turn equals ^\h- So we would have e'^'^" = -^Ih- From Eqs. ()4.61|) and ()4.65p we 
know that /5' > 0, whereas < between origin and horizon. As e^^ equals ^ at 
the origin, the two metric functions cannot attain the same value at the horizon 
unless they are constant functions, which is the case for vanishing coupling 77 = 0. 

Summarizing we get 

= -, e'^- = (4.70) 

One sees immediately, that Ph can only be finite if r/ < 0.5. CSS solutions, 
regular at both origin and past SSH, can therefore only exist for small couphngs 
r] G [0,0.5). 

Note that solutions to ^M>-^!E^ satisfying piTTjl and IKT^ (if they exist) 
are analytic in z and the shooting parameters. (In the vicinity of the singular 
points z = and z^ one can use Prop. 1 of [14^ to show analyticity, for other z 
it follows from the analyticity of the right hand sides of the equations.) 

In order to construct regular solutions numerically, we proceed as follows. We 
consider the boundary value problem, consisting of the coupled system of four 
first order ODEs Eqs. ()4.fjl|l . ()4.f)4|l and ()4.fj5|l subject to the boundary conditions 
fj4.67p and (pn = f • The relation for Ph can be dropped here, as it results from a 
subsidiary equation, which is not used for the calculation. (An alternative would 
be to substitute Eq. ()4.65|) by the algebraic relation ()4.66|) . thereby reducing the 
system to three first order ODEs.) The four free shooting parameters at the 
boundaries are 

0'(O) = 6, ZH, (P'h, Ph. (4.71) 

For fixed 77 the problem is solved numerically using the shooting and matching 
routine d02agf from the NAG Library [^0]. We start at 77 = by pinning down 
one member of the discrete one parameter family of solutions reported by Bizon 
|Hj and follow this solution to higher values of the coupling constant. As we are 
mainly interested in the ground state and first excited state, we only constructed 
the first few solutions. 



4.2.2 Phenomenology of Numerically Constructed Solu- 
tions 

As in [11^ we find that the solutions present on Minkowski background stay 
regular between origin and past SSH up to a maximal value of the coupling 
constant r]max = 0.5. At this maximal value the metric function P diverges at 
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the horizon, which can be inferred from ()4.7()|1 . A more complete picture of this 
phenomenon can be obtained if one considers the analytic continuation beyond 
the past SSH, which will be discussed in the next section. 

Figure shows the shooting parameter 0'(O) for the ground state and the first 
excited state as functions of the coupling. This parameter alone determines the 
solution. Figure shows the value of (3h for the two solutions again as functions 
of the coupling. 



4.2.3 Analytic Continuation beyond the past SSH 

In order to get a more global picture of the solutions described above, we study 
the analytic continuation of the solutions beyond the past SSH. We do this by 
integrating Eqs. ()4.61|) . ()4.64|) and ()4.65|) towards larger values of z, with initial 
conditions imposed at the horizon zh, that are taken from the solutions to the 
boundary value problem discussed above. Note, that - as the past SSH's domain 
of dependence is zero - in general it is not enough to give data there. This 
procedure only works by the means of analyticity. 

As already mentioned above, we have (— )' = at the horizon. From this it 
follows, that 2z — ^ > for z close to but larger than zjj- Therefore (^)' > 
for z > Zh and ^ is monotonically increasing in the region beyond the horizon. 
So it might happen, that ^ eventually equals 2z at some value zs, which means 
that the equations have a second (or counting the origin a third) singular point 
there. 

It was already shown in JT] that this second singular point does not correspond 
to a second self-similarity horizon (i.e. a null surface that is mapped to itself by 
the diffeomorphisms) , but rather a spacelike hypersurface where 2m /r —>■ 1. We 
will repeat the argument here. 

Consider the expression 

h{z) = ^—y. (4.72) 

r 

This function is positive for zh < z < zs and diverges, when approaching the 
past SSH from above: 

lim h{z) = +00. (4.73) 
Using Eqs. ()4.61|1 and ()4.65|1 we find for the derivative 

^'W = (^^TTp (-2(2^ - 7)^' - (2 - (7)')) 



2e'^^ 



(4.74) 
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Figure 4.2: The shooting parameter 0'(O) for the ground state (dashed hne) and 
the first excitation (sohd hne) as a function of the couphng constant rj. 



2.5 




V 



Figure 4.3: The shooting parameter Ph for the ground state (sohd hne) and the 
first excitation (dashed hne) as a function of the couphng constant i]. (3h rises for 
both solutions as one approaches the maximal coupling rjmax = 0.5 and diverges 
in the limit, as can be inferred from Eq. I4.7UI 
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which is less than zero in the region we are interested in. From this it follows, 
that 

lim h{z) exists, (4-75) 

and therefore as the denominator vanishes, we have 

lim_ e-2^(") = 0. (4.76) 

From Eq. ()4.66jl one finds, that the combination (-^ — 2z) e^^^(0')^ stays finite 
in the limit z —>■ zs- Making use of this one can show that the scalar curvature, 
the Kretschmann invariant, the square of the Ricci tensor R^^R^^'^ and the Weyl 
invariant are bounded when z — > zs- Up to now it is not clear whether itself 
has a limit. If so, the above invariants have a limit as well. 

As in [TT] our numerical integration gives the following results: for each member of 
the one parameter family, there exists a critical value of the coupling constant ?7* 
such that for smaller values of the coupling, the solution extends smoothly up to 
the future SSH, whereas for stronger couplings the geometry contains marginally 
trapped surfaces aX z = Zs > z^. The coordinate location of these marginally 
trapped surfaces is a decreasing function of the coupling constant and eventually 
merges with the location of the past SSH in the limit t] — > rjmax- Furthermore 
the critical value of the coupling increases with the excitation number n. 



4.3 Stability of CSS Solutions 

In order to answer the question, whether any of the above constructed solutions 
may play a role as a critical solution in gravitational collapse, it is essential to 
study the stability properties. 

In jH] the stability of the one-parameter family of CSS solutions on Minkowski 
background was studied. The results reported are, that the n-th excitation has n 
unstable modes, in particular the ground state is stable and the first excitation 
has one unstable mode. This suggests, that the ground state plays the role of 
a global attractor in the time evolution of strong enough initial data, and that 
the first excitation might be a critical solution at the border of two different end 
states. These predictions were verified in ^H]- We will talk about these critical 
phenomena in more detail in Sec. 15.71 

As the limit r/ ^ is a regular limit concerning the existence of CSS solutions, 
one expects, that the stability properties of the solutions do not change, when 
gravity is switched on, at least as long as the coupling constant is small. 

We report here on the stability analysis for the coupled case, which we performed 
in two essentially different ways, and show, that indeed the stabihty properties 
of the CSS solutions do not change for small couplings. 
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Figure 4.4: The field of tlie first excitation for r\ = 0.1. Tlie vertical line marks 
the horizon at zh — 0.45457. 
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Figure 4.5: The same situation as above. Plotted are the metric function ^ as 
well as — for the first excitation at 77 = 0.1. — crosses the line 2z at the location 
of the horizon. As r) is well below the critical value rjl ~ 0.152, — stays well 
below 2z outside the horizon and ^ is far from being unity anywhere. 
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Figure 4.6: The field of the first excitation for 77 = 0.175. The vertical line 
marks the horizon at = 0.4198. The right boundary of the plot is at the 
second singular point 2;^ = 0.9292. 
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Figure 4.7: The same situation as above. Plotted are the metric function ^ as 
well as — for the first excitation at 77 = 0.175. — crosses the line 2z at the 
location of the horizon. As rj is above the critical value r)l ~ 0.152, equals 2z 
for a second time at zs — 0.9292. At the same time — tends to 1. 
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Figure 4.8: The coordinate locations of the past self-similarity horizon zh and the 
second singular point zs for the first CSS excitation as functions of the coupling 
constant rj. For rj < rjl c:^ 0.152 the analytic extension of the CSS solution is 
regular up to the future SSH {z — >• cxo). For bigger rj the solution develops an 
apparent horizon at the spacelike hypersurface z = zs = const. This hypersurface 
finally merges with the past SSH in the limit rj 0.5. The other excitations 
behave in a similar way, where ?7* is an increasing function of the excitation 
number n. 

The first method is the standard way to determine the stability properties of 
CSS solutions. The CSS solution is perturbed in a small spherically symmetric, 
time dependent way. The perturbations then are decomposed into modes with 
an exponential (and oscillatory) time dependence. Inserting this ansatz into the 
linearized field equations gives a coupled system of ODEs, with the same singular 
points as for the background solution. Regularity then requires the perturbations 
to be solutions to a boundary value problem, or to be more precise to a linear 
eigenvalue problem, the real part of the eigenvalue, if bigger than zero, being 
responsible for the exponential growth of the unstable mode. The advantage of 
this method is that one only deals with ODEs, which can be integrated rather 
accurately. One main disadvantage is, that if using a shooting and matching 
method, one needs good initial guesses for the shooting parameters. In other 
words, one can never be sure, that one obtains all of the relevant, i. e. unstable, 
modes, unless one has further theoretical arguments. Although in all similar 
situations, where the stability of an "expected-to-be" critical CSS solution was 
analyzed, the unstable modes were all real, there is no theory guaranteeing this. 
For the system on fixed Minkowski background it was shown in ^j, that the 
perturbation operator can be brought into a self-adjoint form (using orthogonal 
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coordinates), therefore the eigenvalues have to be real. Furthermore a theorem 
for Sturm-Liouville operators could be applied, which determined the number of 
eigenvalues giving rise to unstable modes. ^ Unfortunately a similar analysis does 
not exist for the coupled system. 

We therefore present an alternative method to compute the unstable modes. This 
method was proposed and carried out by J. Thornburg [Slj. The method uses 
the full nonlinear field equations. It is based on the observation, that a numerical 
time evolution maps the discretized field, i.e. the N-dimensional vector, on the 
initial slice to a discretized field on a later slice. If the initial configuration is close 
to a CSS solution and the time step is small enough^, then the relation between 
the deviations from the CSS solution at the initial and final slices is linear, i.e. 
determined by an x matrix. The unstable modes can then be extracted 
from the eigenvalues and eigenvectors of this matrix. The main advantage of 
this method is, that it should give all the unstable modes of the CSS solution 
that are "seen" by the time evolution code, and therefore, if the number of grid 
points is big enough, all the unstable modes of the continuum problem. A 
further advantage of this method is, that it uses an already existing evolution 
code (the DICE code, see App. IH}. One minor disadvantage of the method is, 
that the numerical answers are expected to be less accurate than the answers 
obtained from the ODE boundary value problem. Unfortunately up to now this 
method suffers from a more serious drawback, namely it fails to converge with 
respect to resolution. More precisely, increasing the number of grid points, the 
numerical results move further and further away from the predicted ones (an 
example is the gauge mode, described below). What is even worse, is that we 
do not know, why this method does not converge. One possible reason could be 
that for higher resolutions the method gets increasingly ill-conditioned, that is 
increasingly sensitive to small numerical errors in the numerical time evolution. 
Nevertheless, as we get answers for a certain number of grid points (A^ = 500) 
which are in good agreement with the results of the other method, we are inclined 
to believe these results (for this number of grid points) and deduce the stability 
properties from them. 

^We mention that the continuous spectrum of the operator in |S] would seem to be unstable 
in adapted coordinates (t, z). In fact these modes are growing as fast as the gauge mode, 
which shows, that the growth is due to the shrinking of adapted coordinates and therefore 
these modes are not considered as unstable. We also mention that members of the continuous 
spectrum oscillate infinitley many times in the vicinity of the horizon. They therefore cannot 
be detected by the methods described below. 

^The time step under consideration has to be small such that unstable modes don't drive 
the solution out of the linear regime. The time step may consist of several numerical time steps 
of the evolution code. 
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4.3.1 Unstable Modes from a Boundary Value Problem 



In order to analyze the stability of the CSS solutions, we proceed in the usual 
way. Consider small time-dependent radial perturbations of the CSS solution 



0(r,z) = (j)n{z) + S(f){T,z), 
j{r,z) = {jUz)+6V{t,z), 



(4.77) 
(4.78) 

(4.79) 



where (pn, Pn, i^)n denote the n-th CSS excitation. (Note that by 5V we denote 
the perturbation of the metric function — not of V alone). As the perturbations 
are supposed to be small, we can linearize Eqs. ()2.48|) . ()2.55|) . ()2.56|) in these 
perturbations. Together with the requirement that the background solutions 
(pn, Pn, {^)n solvc the CSS equatious ()4.(il|l . ()4.(j4j) and ()4.()2j) (or alternatively 
fl4.(j5j) ). one obtains the following linear system of PDEs: 

6(3' {t,z) = zr]^'^{z)64>'{T,z), (4.80) 

6V'iT,z) = -^(2e2'3"W (-l + 2r/sin(0„(z))') (5/3(r,2) + (5V(r,z) + 

+ 2e^'^"^'^r] sm{2(f)n{z))6(j){T,z)], (4.81) 



2e2/3„(^) l2z-i-Uz)) 6(3iT,z) sin(2 </)„(z)) + ^ (-1 + 2 sin(0„(z))') 0;(^) 



+ 5V{t,z) [2z2<^'„(z) + e2/3"W (sin(2 0„(z)) + ^ (-1 + 2r7sin(0„(z))') 0;(^)) 



+ 



+ [2z-{-)4z) (2e^('"^'U4>{r,z) [cos{2 ^z)) + z v sm{2 Mz)) 4>'niz)] + 



+ z 



2\ V 

sin(0„(z)) ) - {—)n{z) 



50'(r,^) + 



+ z[2z- {-)n{z) ) 50" (r, z) + 2 (50(r, z)\] + 



+ 



I 22-(-)„(z) 
r 



y{r,z). 



We now decompose the general perturbations ()4.77|) into modes of the form 



e^^y{z) 



5(I){t,z) : 
6P{t,z) : 

SViT,z) = e^^h{z). 



e'-^giz), 



(4.83) 



Inserting this ansatz into Eqs. ()4.8()j) . ()4.82j) yields a coupled system of ODEs 
for the perturbations y{z),g{z) and h{z). This system again suffers from two 
singular points, the origin z = and the past SSH z = zh- 



(4. 
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Regularity at the origin restricts the perturbations to 

y(0)=0, 9(0) = 0, h{0) = l, (4.84) 

whereas the first spatial derivative of the field y'{0) is unconstrained. 

At the horizon the requirement of regularity relates y'^j to the boundary values 
of the other fields 

= X2zh{-1 + 2r/) (^~^ + ~ ^"^^^^^^ ^ 4zh{<P'Jh{-1 + 2v)gH - 
- i<P'jHil + 2v + zlr^il-2r^)i<P'X)hH). (4.85) 

Eqs. (HSni), (fO^ together with the boundary conditions P^Mjl . consti- 
tute an eigenvalue problem, with eigenvalue A and eigenfunctions {y,y',g, h). As 
already mentioned the eigenvalues and eigenfunctions may be complex. Never- 
theless, as the coefficients in Eqs. ()4.8U|) . ()4.82|) are all real the eigenvalues and 
eigenvectors, if complex, come in complex conjugate pairs. 

The resulting problem again is a boundary value problem, consisting of eight 
linear first order ODEs (Eqs. ()4.80| ()4.82|l ) separated into real and imaginary 
part) and the eight boundary conditions (real and imaginary parts of Eqs. ()4.84|1 
and ()4.85j) ). The parameters that have to be matched are 

A, yH, Qh, hu, (4.86) 

(again real and imaginary parts thereof). As the problem is linear and homoge- 
neous the solutions are only fixed up to an overall scale. We fix that, by setting 
l/'(0) = 1. 

This boundary value problem again was solved numerically, using the shooting 
and matching routine d02agf of the NAG-library ([HD])- The background solu- 
tion was computed first and interpolated when necessary in order to provide the 
coefficients in Eqs. fOn|l . (g^l. 



4.3.2 Gauge Modes 

Before reporting on the numerical results, we determine the gauge modes, that 
result from a certain arbitrariness in relating the adapted coordinates (r, z) to 
Bondi coordinates (M,r). Recalling from Sec. 14.1.21 (Eq. ()4.39|l ) we have 

t{u) = — ln(u* — u) u{t) = u* — e^^ 

J. ^ 

^iu,r) = — r{T,z) = e ^z. (4.87) 

[u* — u) 
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A shift in the culmination time u* u* = u* + e (corresponding to a shift of the 
origin of u) yields another pair of adapted coordinates {f,z), defined as above 
with respect to u*. The relation between these two sets of adapted coordinates 
is given by 

f = 
z = 

For small e we linearize to get 

f = T — ee'^ = T + ex'^ {t, z) 

z = z - eze" = z + ex\r,z), (4.89) 

where we have introduced the generating vector field x'^(r, 2;). In the following 
we will treat the vector field ^is general and show in the end, that the form 
()4.89j) is indeed the only possible one.^ 

Such a small coordinate transformation introduces a (small) change in the per- 
turbations according to 

Stj) = 5<p - eCy,(t)Q, (4.90) 

where in our case the objects with the subscript denote a CSS solution. Pure 
gauge modes are characterized by 

5<P = eC^(t>o. (4.91) 

i.e. modes of this form can be removed by a small change in the coordinates 
according to ()4.89j) . 

We first note, that the coordinate transformations are not completely arbitrary, 
but have to ensure that a hypersurface f = const still is a null cone, which is 
reflected by 

^x(^?o).. = 0, (4.92) 

or 

X:, = ^x^ = x^{r). (4.93) 

The second observation is, that via 14. 771 we fixed the "S^ - part" of the metric to 
be e~^'^2;^, which gives 

CM)ee = 0, (4.94) 

*In particular the following arguments rule out a coordinate transformation, that resets the 
origin of r, r — > r + c : such a coordinate transformation violates (|4.98() and therefore does 
not correspond to a regular perturbation of the CSS solution within our choice of coordinates 
(especially the choice of u being proper time at the origin) . 



In 1 



ee 



1 + ee^ 



(4. 
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or 



x''{9o)ee,T + x'{9o)ee,z = 0, 



(4.95) 



For the remaining components of ()4.9H) we get 



-,-2t 



2{-x^ + x:,)e2*((^)o - 2z) + r{e'^-{{^)o - 2z))' + 



2(-x^ + x:.)e2*((^)o - 2z) + zx\e^^\{^-\ - 2z))' + 



+2zx]re 
rU^) = zx^^'oiz). 



(4.96) 



On the other hand the perturbations of the metric functions P and — are re- 
lated to the perturbations of the components of the metric with respect to (r, z) 
coordinates via 



.e-2rg2/3o 



2{{-)o-2z)5p + SV 



-2e-2^e2'^°<5/3. 



(4.97) 



The regularity requirements at the origin discussed in the last section, namely 
5f3{T,0) = 6V{t,0) = 0, give ^(^^-^('r, 0) = (5(7t-z(t, 0) = 0. Comparing this to 
()4.9fij) gives the further restriction on the generating vector field 



(4.98) 
(4.99) 



and therefore 

= e^ x" = ze^, 

which is exactly the coordinate transformation introduced above ()4.89|) . 

Finally we combine ()4.96|) . ()4.97p with ()4.99|) . In a stability analysis, as described 
above, we therefore expect to find the gauge mode 



5p 
6V 



ee^z(J)L 



ee z — L 



(4.100) 
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4.3.3 Unstable Modes from a Matrix Analysis 



The following method was proposed and carried out by J. Thornburg [TU]. It uses 
the full (nonlinear) time evolution equations. As already mentioned, the method 
has the advantage, to give all the unstable modes, in contrast to the shooting 
and matching method described above, where only those unstable modes can be 
found, which lie close to the "initial guess" . 

Consider Eqs. (j23Hl), and ^H^ . As described in Sec. O it is sufficient 

to prescribe the matter field at the initial null cone u = uq (or r = tq). The 
metric functions on the initial slice then are determined by via the hypersurface 
equations ()2.55p and ()2.56|) . Eq. ()2.48|) is then used to evolve the configuration 
(together with Eqs. ()2.55|1 and ()2.56|) in order to update the metric functions). 
Altogether the time evolution maps the field 0(to, z) at the initial slice to a field 
configuration 0(to + At, z) at a later slice: 

0(ro + Ar,^) = F(0(ro,^)), (4.101) 

where F is the nonlinear operator representing time evolution. Of course the 
operator F depends on the time step Ar. Clearly, if the initial data correspond 
to the CSS solution, F acts as the identity operator, 

F{(t)css{ro, z)) = 0c55(r, z). (4.102) 



In an actual numerical time evolution in (1+1) dimensions any smooth function 
of z is represented by it's values at the grid points, i.e. the function (j){TQ,z) is 
replaced by the N-dimensional vector (0°)j=i^7v, where 0° = 0(ro,2;*), and is 
the number of grid points. The numerical time evolution then maps this vector 
to the corresponding vector at the next time step: 

= F,{ct>]). (4.103) 



Consider now a small generic perturbation of the CSS solution, (pcssijf)) z) + 
50(ro,-2). If the time step Ar is small enough, this configuration is mapped to 
another small perturbation of the CSS solution at time tq + Ar. Translating this 
to the language of a finite grid we have 



0! = (0CS5)! + (50)! = i^.((0CS5); + ('50);). 

For a small perturbation this can be linearized to give 



(0C55).' + (50)!=i^.((0 



css)j) + 



dF, 



m' = (0C55)° + 



OF 



(•Piss) 



I'CSS' 



(4.104) 



(4.105) 
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So for the perturbation the following linear relation holds 



fcss 



{6<P)'j. (4.106) 

) 



The Jacobian on the right hand side is a x matrix, which depends on the 
size of the time step At. Therefore it has eigenvalues {Xi}i=i^N, which are 
functions of the time step At. If the Jacobian is diagonalizable we can switch to 
a basis such that 

= Ai((50)°, with Ai = A,(Ar). (4.107) 
For the quotient of differences we get 

6(j)i{To + AT) -6(j),{To) Ai-1 fA,no\ 

^ = -^^Mro). (4.108) 

In order that the limit At exists, the denominator on the right hand side 
has to be proportional to At. We set 

Xi-1 = \i At, (4.109) 

Aj being constants. 

In the limit At ^ we have 



and therefore 



,(r) = A,50,(r) (4.110) 



r) = e^^(^-^°)5(/.,(ro). (4.111) 



This is valid with respect to the eigenbasis of the Jacobian. Transforming back 
to the original basis we have 

{Scf^k)^{r) = e"'^-(^-^«)(50,),(ro), (4.112) 

where Scpk denotes the k-th eigenvector of the Jacobian. The perturbation ()4.112|1 
now is of the same form as the modes in Eq. ()4.8Hj) . with the eigenvalues A being 
related to the eigenvalues of the time evolution Jacobian via Eq. ()4.109|) . 

The Jacobian was computed using the DICE code (described in Appendix Ej). As 
this code uses grid points, that are freely falling along ingoing null geodesies, the 
evolved field had to be interpolated to give the appropriate values at constant 
z. Technically the Jacobian is computed, by first evolving the CSS solution for 
one time step, and then perturbing each grid point separately, according to the N 
perturbations (50^)° = e6ik, k = 1,N, with e small. Each of these perturbations 
is evolved again for one time step, so in the end the N x N numbers Fi{(j)k) are 
known. The Jacobian results from this by a forward differencing. In order to 
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compute the eigenvalues and eigenvectors of the matrix a hnear algebra package 
(EISPACK [67J) is used. 

Finally we want to point out, that although this method only yields modes of 
the perturbation operator, according to the N grid points involved, we expect to 
find all the relevant unstable modes. The reason for this is, that we expect the 
unstable modes to vary on scales that are large compared to the grid spacing. 

4.3.4 Numerical Results of the Stability Analysis 

Before we present our results here, we want to stress two points: First we are 
interested in the stability of the CSS ground state and the first CSS excitation, 
because we want to get information on the possible roles they might play in 
the context of critical phenomena. With respect to this, only the stability for 
small couplings (77 < r7*(0) ~ 0.69 for the ground state and r] < 0.2 for the first 
excitation) is of interest. 

Second we trust the results of the boundary value problem Sec. 14.3.11 On the 
other hand, as already mentioned, the numerical scheme of the matrix analysis 
Sec. 14. 3. 31 does not show convergence with resolution, and results should be taken 
with some care. Nevertheless, as we will show below, for = 500 grid points 
(and T] not too large), the results for the gauge modes and the unstable mode for 
the first excitation are in very good agreement with the theoretical predictions 
and the results of the boundary value problem. Therefore - for these couplings 

- we are inclined to trust the results of the matrix analysis, in particular the 
number of unstable modes. 

Given these caveats, we numerically find that the ground state is stable, whereas 
the first excitation has one unstable mode (both results for 77 not too large (See 

Figs, mm and mni)). 

Fig. 14.91 compares the gauge modes - obtained via the boundary value problem 
of Sec. 14.3.11 - of the ground state and first excitation at a coupling 77 = 0.1 to the 
theoretical predictions y gauge = ^0n- Note that in this figure and in the following 
ones the overall scale is chosen arbitrarily. 

Fig. 14.101 compares the eigenfunctions of the gauge mode and the unstable mode 

- obtained via the matrix method of Sec. 14.3.31 - of the first CSS excitation for 
?7 = 0.1 to the results of the boundary value problem. The agreement is very 
good. 

Figs. 14. Ill and 14. 121 show the real parts of the first few eigenvalues of the ground- 
state and first excitation. The agreement of the results of the matrix method and 
the boundary value problem is good for the groundstate at low couplings and 
very good for the first excitation. 

Finally, note that the inverse eigenvalue of the unstable mode of the first CSS 
excitation 1/Ai as a function of coupling rj is very well approximated by a straight 
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Figure 4.9: The gauge modes for the ground state and the first excitation at a 
couphng ?7 = 0.1. Plotted are the predicted functions, y gauge = zcp'^, (sohd hues), 
and the corresponding results of the shooting and matching method ("x", "+"). 
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Figure 4.10: The eigenfunctions of the unstable mode (A = Ai = 5.58463) and the 
gauge mode (A = 1) for the first CSS excitation at a coupling r] = 0.1. Plotted are 
the result of a 500 points matrix analysis (dots, only every 5th point is plotted). 
These are compared to the results from the boundary value problem (lines): the 
gauge mode is compared to the predicted eigenfunction ygauge{z) = z(t)[{z) and 
the unstable mode is compared to the result y{z) of the boundary value problem 
described in Sec. 14.3.11 
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Figure 4.11: Dots represent the real parts of the first few eigenvalues of pertur- 
bations of the groundstate as obtained by a matrix analysis with = 500. The 
line A = 1 represents the predicted eigenvalue of the gauge mode. For small 
7] the matrix analysis gives no unstable mode. For bigger r] (e.g. at rj = 0.4) 
the additional positive eigenvalues could not be confirmed by the shooting and 
matching method. 



Re{Xi 




Figure 4.12: The same situation as in Fig. 14.11) for the first CSS excitation The 
dashed line starting at r] = with 6.33 represents the eigenvalue of the unstable 
mode obtained from the boundary value problem. The agreement of results of 
the two methods is very good. The matrix analysis gives one unstable mode for 
Tj < 0.35. The additional positive eigenvalues for larger r] (e.g. r] = 0.4) could not 
be confirmed by the shooting and matching method. 
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l/Ai(T/),/(r7) 




Figure 4.13: Plotted is 1/Ai (dots), where Ai is the eigenvalue of the unstable 
mode of the first CSS excitation, obtained by the shooting and matching method 
of Sec. 14.3. H as a function of rj. This is very well fitted by the straight line 
/('?) = 0.219977] + 0.15736 (solid line). See the next figure for the error of this 
fit. 
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line (see Figs. I4.1HI and 14.141 



4.4 Numerical Construction of DSS Solutions 

This section deals with the numerical construction of DSS solutions to the self 
gravitating nonlinear o model. Due to the (periodic) time dependence of the 
metric functions and the field, the numerical construction of DSS solutions is 
considerably more involved than the corresponding construction of CSS solutions. 
Concerning self-gravitating matter fields there are essentially two papers (not 
counting subsequent ones using the same methods) that deal with the problem 
of constructing time-periodic solutions to a boundary value problem (in space). 
Seidel and Suen construct solutions to the Einstein-Klein-Gordon system with 
mass, which are oscillating (periodic in time). Gundlach jHSj constructs a DSS 
solution to the massless Einstein-Klein-Gordon system. Both methods use Fourier 
series, as suggested by the periodic time dependence, but their "implementations" 
are different, as will be explained below. 

Following Gundlach closely, we present a method, that involves discrete Fourier 
transform, pseudo spectral methods and the reduction to a boundary value prob- 
lem for ODEs. 

The procedure of construction can be roughly outlined as follows: we are looking 
for solutions to Eqs. ()2.48|1 . ()2.55j) and ()2.5(i|l . rewritten in adapted coordinates 
fj4.38|) . that are regular between origin and past SSH and that are periodic in 
the adapted time coordinate r. We expand the metric functions as well as the 
field into Fourier series in time, where the coefficients are functions of the spatial 
variable z. Inserting these Fourier series into the equations yields a coupled sys- 
tem of ODEs for the Fourier coefficients. This system has singular points at the 
origin and the past SSH, so together with the boundary conditions required by 
regularity we have to solve an ODE boundary value problem for the Fourier co- 
efficients. As the Fourier series consist of infinite many terms and as the problem 
is nonlinear, one has to truncate the Fourier series at some appropriate maxi- 
mal frequency. This then yields a boundary value problem for a finite system 
of ODEs, which again can be solved by the means of a shooting and matching 
routine. 

Seidel, Suen and Gundlach use different methods for explicitly setting up the 
ODEs. Seidel and Suen truncate the series and plug the truncated series into the 
equations (see Appendix |Bj) . Comparing coefficients yields the desired coupled 
system of ODEs for the coefficients. The disadvantage of this method is, that 
these direct expressions quickly get horribly complicated as one allows for more 
and more coefficients. 

Gundlach on the other hand uses so called pseudo-spectral methods, that is, he 
does part of the computations in "real" space and part in Fourier space (see 
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Appendix El) • Basically the variables are the Fourier coefficients, but in order to 
set up the ODEs, the coefficients are transformed back to real space, there the 
algebraic manipulations are carried out in order to define the derivatives of the 
functions, which then in turn are transformed back to Fourier space. This requires 
a pair of (backward and forward) Fourier transformations at each integration 
step. Employing a fast Fourier transform (FFT) instead of the ordinary discrete 
Fourier transform (DFT) is essential for reducing computational cost. 



4.4.1 The DSS Equations 

We start by transforming Eqs. ()2.48|) . ()2.55|) and ()2.56|) to the adapted coordinates 
defined in Sec. 14.1.21 This gives 

= 1^(0')' (4.113) 



(-)' = —{-e'^ + 2r]e'^sm\(j)) + -] (4.114) 
r z \ r J 

and 

(j)" = 7 ^ (-e^^ {sin(2 0) + z (-1 + 2/5sin(0)') 0'} + 

+ ^ |(^^+4^C-4C^)0'-2C(0 + ^0')}), (4.115) 

where ' = dz and ' = dr- Looking for DSS solutions means that we require 
P = P{t,z),^ = ^{T,z),(f) = (I){t,z)X = CiT)^ where all the time dependencies 
being periodic, with some period A, the determination of which is part of the 
problem. 

Furthermore we fix the coordinate freedom contained in the function (, by 

C(r)-C(r) + ^-(r,l) = 0, (4.116) 
2 r 

which makes the hypersurface z = 1 null, as explained in Sec. 14.1.21 Assuming 
that ^(r, 1) is given, Eq. ()4.116p is a first order ODE for C('r)) with periodic 
boundary conditions. 

As we will have to deal with an equation similar to Eq. ()4.116j) again, we have a 
short look at the more general equation 

f(r)+g{T)f{T) + h{T) = 0, (4.117) 

where g{T) and /i(r) are given functions, periodic in r with period A, and /(r) 
is the unknown, which is required to be periodic too with the same period. 
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The general solution to Eq. ()4.117|1 is given by 

\ TO / T"0 

where the function G{t) as defined above was introduced for abbreviation. The 
constant c entering the formula now has to be determined from the required 
periodicity of /(t). The behavior of G{t) under a shift of A is 

G(r + A) = G(ro + A) + G(r), (4.119) 

so 

To+A 

/(r + A) = /(r) - c e-^W (l - e-^(^°+^)) - e-^(-)e-^(^°+^) / e^(^)/i(f )rff . 



TO 



(4.120) 



And so, if G(ro + A) 7^ 0, the constant is determined to be 

To+A 

e^^^^h{f)df, (4.121) 



1 

~ 1 _ eG(To+A) 

To 

and therefore the unique periodic solution to Eq. ()4.117j] is given by 

/ T ro+A \ 

fir) = e-^(^) - I e''^'^h{f)df + ^ _ J^^^^^^ J e''^'^h{f)df . (4.122) 



In case gi^r) has no zero frequency, i.e. G{tq + A) = 0, there is no periodic solution 
to Eq. ()4.117|1 unless the last integral in ()4.12()j) vanishes as well, in which case 
there is a one-parameter family of periodic solutions to Eq. ()4.117|1 of the form 

According to this the solution to Eq. ()4.11f)|l can be given in closed form: 

(T TQ+A 
- / e-^(r, l)df + ^J-^ J e-^(f , l)df | . (4.123) 
TO TO 

There is an alternative way to compute (, namely as Eq. ()4.11(j|l is linear in all 
periodic functions, the Fourier coefficients of ( can be computed directly from 
the Fourier coefficients of — (r, 1). We will give details on this in Sec. 14.4.41 
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Finally we rewrite Eqs. ()4.113j) - ()4.115j] . using the logarithmic coordinate y = 
In z instead of z. This is useful for numerical purposes, as the main variations in 
the solutions occur close to the origin. We have 

= (4.124) 
(jY = -(-e^^ + 2T]e^f'sm\^) + ^^ (4.125) 



and 



-^ + 2eyc-2eyc 



(-e^^ {sin(2 0) + (-1 + 2 sin(0)^) 0'} + 



+ |(2eH-2Ce^)0'-2Ce2'(0+ 0')}) , (4.126) 
where ' now denotes dy and y = ln(2;). 

Note that r does not enter Eqs. (I4.113jl - (I4.115jl . f (14. 12411 - (14.1 2 (il) respectively) 
explicitly. Therefore, given a solution, abbreviated by Z{t,z), all expressions 
resulting from this by a constant shift in r, Z{t + const, z) are solutions to the 
system as well. 



4.4.2 Regularity at Origin and Past SSH 

Regularity at the origin 2; = 0,?/ = — ooisas usual given by 

P(t, z = 0) = 0, -{t,z = 0) = 1, cj){T, z = 0) = 0, (4.127) 

r 

and </)'(t, 2; = 0) is a free (periodic) function of r. Using the logarithmic coordinate 
y = ln{z), we start our numerical integration at some finite j/q. Near the origin 
z = (i.e. y —00) the functions behave as 

f3{T,y) = 0{e'y) 

-(r,y) = Oie'y) 
r 

ct){r,y) = 0i(r)e^ + O(e2^) 

0'(r,y) = 0i(r)e^ + O(e2J'), (4.128) 
where 0i(r) = dz(t){r, z = Q). 

At the past SSH z = 1 (i.e. y = 0) the denominator in ()4.12fij) vanishes. Regularity 
therefore enforces the nominator to vanish as well, that is 

= -e^^^ {sin(2 0H) + (-l + 2r/sin(0^,)') 0'^} + 

+ {((^)/^ + 4C-4c)0'^-2c(0// + 0'h)|, (4.129) 
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where the expressions with subscript H are periodic functions of r. Suppose now, 
we are given the functions Ph,{^)h and (pn, then Eq. ()4.129|) is an ODE in time 
r for In fact, if we abbreviate 

^ = ^((^)/^ + e^'^-(-l + 2r^sin2 0^)^ (4.130) 

h = ^e^fe sin(20H) + (4.131) 

and use Eq. (j4.116p to replace (, then Eq. ()4.129|) is of the same form as Eq. ()4.117|) . 
described in the last section, and has the unique periodic solution ()4.122|) . 



4.4.3 Additional Symmetry 

From critical searches (see Chap.Elfor details), we know, that the critical solution 
for large couplings, which is DSS, not only is periodic with period A, but has an 
additional symmetry, namely 

PiT + A/2,z) = PiT,z), 

V V 
-(r + A/2,^) = -{t,z), 

ar + A/2,z) = C(r,z), 

0(r + A/2,2) = -0(r,2), (4.132) 

which means that the metric functions as well as ( consist only of even frequen- 
cies^, whereas the field and it's derivatives contain only odd frequencies. As we 
are interested in the direct construction of the DSS solution, which is the criti- 
cal solution in a certain range of coupling constants, we impose this additional 
symmetry on the solution by the means of its Fourier coefficients. 



4.4.4 Numerical Construction of DSS solutions via an ODE 
boundary value problem 

According to the required periodicity in r of the solution we expand the metric 
functions, (, and the field into Fourier series. The truncation of theses series 
is performed according to the discrete Fourier transform described in detail in 
Appendix El We denote the number of "collocation points" in "r space" by N, 
assuming, that it is an integer multiple of 4. So the equally spaced points in r 
space are given by = kA/N. 

In principle these N points give rise to Fourier coefficients, but the additional 
symmetry introduced in the last section. Sec. 14.4.31 reduces this number by a 

^even multiples of 1/A. 
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factor of 2. So each function consists of M = N/2 Fourier coefficients, wliere we 
assume, that M is even. The expansions then are given by 

pN{Tk,z) = po{z)+ ^ {pcos)2i{z) cos{^^) + ^ {(3 sin) 21 (z) sm{^^) + 

1=1 1=1 



+{(3cos)Miz) cos(7r/c) 

M/2-1 ^ M/2-1 ^ 

—)N{Tk,z) = T^o(^)+ 2^ (V^cos)2K^)cos(-^) + 2^ (T/szTijaK^) sm(-^) + 

Z=l i=l 

+ (Vcos)m(-2) cos(7r/c) 

Civ(Tfe,^) = Co(^)+ (Ccos)2K^)cos(-^) + (Cs^^)2K^)sin(-^) + 

1=1 1=1 
+ {(cos)m{z) cos(7rA;) 

(4.133) 



and 



M/2 r, /or iN( -'^/^ 



J. / N V-/^ N / N M2l-l)k. , , , .2TT{2l-l)k. 

(t>N[Tk,z) = 2^(0cos)2i-i(^)cos( — )+2_^{(t)sin)2i-i[z) sm( — ), 

1=1 1=1 

(4.134) 

and the expansion for 0' follows directly from the one of 0. The coefficients for 
a function with only even frequencies are given by e.g. 

N-l 

= nT.P^^^^) (4.135) 



A;=0 



2 ^ ATTlk 

{Pcos)2i = i = h-- ■M/2-1 (4.136) 

fe=0 



(/3szn)2^ = ^ = 1,-- ■M/2-1 (4.137) 

A;=0 



(/3cos)m = — J]/9Ar(Tfe)cos(7rA;), (4.138) 



TV 

A;=0 



and for a function with odd frequencies by 

N-l 



{ct>cos)2i-i = 1^ E 0iv(r,) cos(^^^^^^^) Z = l,...M/2 (4.139) 

k=0 

{(t)sin)2i-i = ^5^0^(rfc)sin(^^^^^^^il^) / = !,... M/2. (4.140) 

fe=0 
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Our variables that have to solve the coupled system of ODEs now are the AM 
Fourier coefficients of /3, ^, and 0'. Note that ( is not part of the system, as it 
can be computed whenever — is given at the horizon. The boundary conditions 
consist of boundary conditions for the coefficients of /3, ^ and (f) at the origin, 
and on the other hand of the conditions on the coefficients of 0' at the horizon, 
so in total we have 4M boundary conditions imposed on our system of 4M ffist 
order ODEs. 

The free (shooting) parameters consist of the coefficients of /3, ^ and at the 
horizon, which make a total of 3M. Furthermore the M coefficients of 0i(r) = 
dz4>{T,z = 0) are free parameters. Nevertheless, as noted in Sec. 14. 4. Ti the equa- 
tions ()4.113|1 - ()4.115|) are translation invariant in r. A constant shift in r therefore 
transforms a given solution again to a solution. On the other hand, a constant 
shift in r just changes the Fourier coefficients in a well defined way. We can 
therefore chose one coefficient of 0i (r) arbitrarily, so we are left with only M — 1 
shooting parameters at the origin. Finally there is the period A, which is the 
last free parameter. So in total we have again 4M shooting parameters. 

In order to solve the boundary value problem for the ODEs again a shooting and 
matching method is used (see Appendix El- We describe now in detail, how the 
system of ODEs is "set up" numerically. The ffist step consists of providing good 
initial guesses for the shooting parameters. Given these guesses, Eq. ()4.116|) is 
solved for (. As this equation is linear in the periodic functions, it can be solved 
directly in Fourier space. Applying the rules for differentiation in Fourier space, 
as given in Appendix El to and comparing coefficients, Eq. ()4.11(jj) has the 
following solution in Fourier space 

Co = ^m,=o (4-141) 

1 / Airl \ 

iCs^nU = 2(l + (4vr//A)^) (^^osU^^o+ {Vs^nU^^,j4.U3) 

{Ccos)n/2 = liVcos)N/2\y^o. (4.144) 

In the next step, the coefficients of (p'^j have to be computed from the coefficients of 
( and the other variables at the horizon. To do this, the variables are transformed 
back to real space. There formula ()4.122j) is used to obtain 0'^^ as a function of r. 
Transformation to Fourier space then yields the desired coefficients. Note, that 
each time such a "backward-forward transformation pair" is used, the number 
of coefficients is ffist increased by some factor and after the operations in real 
space have been performed and the forward transformation has been applied, 
the number of coefficients is reduced to the original size again. This is a way to 
reduce aliasing errors, as explained in more detail in App. El 
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After these first steps are completed, all the variables, i.e. the Fourier coefficients, 
are known at the boundaries, the coefficients of ( are determined and the period A 
has some definite value. In order to integrate the ODEs, the spatial derivatives 
of the Fourier coefficients have to be determined. This is achieved again by a 
transformation to real space. There the right hand sides of Eqs. ()4.124j) - ()4.12(i|l 
are evaluated, and a transformation back to Fourier space yields the desired 
derivatives of the coefficients. Again for de-aliasing the number of coefficients is 
first increased and after all the operations reduced again. These operations have 
to be carried out at each integration step, which is determined internally by the 
NAG routine's integrator. This huge amount of Fourier transformations for a 
single integration to the matching point, necessitates the use of the fast Fourier 
transform. 

The shooting and matching routine d02agf again was taken from the NAG li- 
brary (pU]). We mention one technical detail, concerning the magnitude of the 
variables. As the solution is expected to be smooth, the "higher" Fourier coef- 
ficients should decrease exponentially in magnitude. As the NAG routine uses 
some mixture of absolute and relative error to determine the local error of the 
solution, it is necessary to rescale the variables to approximate unity. 

As the shooting and matching method needs a good initial guess, we use the 
results of a "dice-critical-search". (For details on the setup of a critical search 
and on the following see Chap. El). As our initial guess we take the critical 
solution at ?7 = 100, which is clearly a DSS solution. Since these data are given 
in terms of Bondi coordinates we have to switch to adapted coordinates. We 
compute u* and the initial guess for the echoing period A from max(2m/r). The 
past self-similarity horizon theoretically is determined as the backwards light 
cone of the culmination point. Another way to determine the SSH is to look for 
the ingoing null geodesic, along which the metric functions /3 and ^ as well as 
the field are periodic functions of r. This gives us the initial guesses for the 
shooting parameters at the horizon. From ^ at the horizon, we compute ( and 
finally dr(f){T,r = 0) can be converted into 0i(t), as the relation between Bondi 
coordinates and adapted coordinates is already fixed. 

We start at r/ = 100 with = 16, and follow the solution down to smaller values 
of 1]. As will be reported in the next section, the echoing period A increases 
"sharply" for rj below 0.5. This increase of the period will necessitate a larger 
and larger number of coefficients. 

Given a solution obtained with a certain number of Fourier coefficients M, we 
can double the number of coefficients with the second half set to zero. This way 
we obtain a reasonable initial guess for the problem with 2M coefficients. 
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4.4.5 Numerical Results 



Applying the method described in the last section, Sec. I4.4.4[ we find good nu- 
merical evidence, that the system Eqs. ()4.124|) - ()4.126p admits a discretely self 
similar solution for 100 > rj > 0.17. The smallest value of t], where we con- 
structed a DSS solution, was r] = 0.17262. At this coupling constant already a 
large number of coefficients is necessary. Lowering t] further would require at 
least = 256, i.e. 128 Fourier coefficients per dependent variable. With this 
number of coefficients - in addition for de-aliasing the number of coefficients was 
increased and decreased by a factor of 8 - a single Newton iteration on an Alpha 
(ev6 processor) already takes several days. 

We find the following behavior: Fig. 14.151 shows that below rj ~ 0.3 the period 
A(?7) rises sharply with decreasing rj. In Sec. 14.5.11 we will give an argument 
that we expect A{r]) to behave like — aln(?7 — rjc) + b for rj close to a critical 
coupling constant rjc- Fig. 14.161 shows a fit of A(?7) against this function for 
0.1726 < 1] < 0.195. According to this fit the period A(?7) would blow up at 
7] ~ 0.17. 

The rise in A (77) also means, that more and more coefficients are needed in order 
to represent the solution to a given accuracy (see Sec. I4.4.6|) . Figs. 14. ITl and HTSl 
illustrate this fact. Fig. 14.171 shows that at a fixed coupling the coefficients decay 
exponentially with the coefficient number. Nevertheless the slope of the decay 
decreases with decreasing coupling. Fig. 14.181 shows the coefficients of the field 
as functions of the coupling rj. With decreasing rj the coefficients grow. 

4.4.6 Convergence with the Number of Fourier Coeffi- 
cients 

As a test for accuracy, we can use the "supplementary" combination of the Ein- 
stein equations, Eq. ()2.57|) . which for an exact solution should evaluate to zero. 
Transformation to self-similar coordinates and multiplication by e~'^({T) yields 



V . V 



2/3 _ _ 9^^2/3 ■2(<^ , 



= 2(3-C C-(C-C)e'^ 2r]e"'sm'((f)) + r]- 

ly ijn ly ry 

- 4r/e^C00'-2r7e^(C-C)(0')^l + 2 0^ - 2r/- 

J r 



(4.145) 



where ' = dy. 



To compute the right hand side of ()4.145|) . the Fourier coefficients of the numer- 
ically computed solution and their derivatives with respect to r are taken, and 



(I4.145|l is evaluated in r space. A Fourier transformation yields the coefficients 
of the expression. Again, in order to diminish aliasing errors, the coefficients 
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Figure 4.15: The period A of the DSS solution as a function of the coupHng 
constant rj. The number of coefficients used to produce the results of this figure 

were: M = 8 (A^ = 16) for 10 < < 100, M = 16 {N = 64) for 0.2933 < r/ < 10, 
M = 32 (iV = 64) for 0.2079 < 77 < 0.2933 and M = 64 (AT = 128) for 
0.1726 <ri< 0.2079. 
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Figure 4.16: "+" represent the period A of the DSS solution in the region of 
smallest rj, where we constructed the solutions. All the solutions are obtained 
with 64 coefficients (A^ = 128). In the region 0.1726 < rj < 0.195 these data were 
fitted against the function f{ri) = —aln{ri—ric) + b. The fit determines the critical 
couphng to be rjc — 0.17, and the constant a — 0.36278. The fitted function f(rj) 
is plotted as dashed line. On the right axis the error is plotted. Presumably due 
to an insufficient number of coefficients the error increases towards the lower end 
of the interval. 
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Figure 4.17: This figure illustrates the exponential fall-off of the Fourier coeffi- 
cients. Plotted is the maximum over y = \n{z) of the Fourier coefficients of 
versus the coefficient number k for the two coupling constants r] = 0.2933 and 
7] = 0.1854. The solutions were computed using M = 32 Fourier coefficients. For 
both couplings the magnitude of the coefficients decreases exponentially. How- 
ever, the slope of this decrease is steeper for larger couplings. 
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Figure 4.18: The magnitude of the Fourier coefficients as a function of the cou- 
pling constant rj. All solutions were computed using M = 32 Fourier coefficients. 
Plotted is again the maximum over the spatial variable y = \n{z) of the coeffi- 
cients of 0. See also Fig. 14.251 for a comparison of errors at the lower and upper 
end of the rj interval shown in this figure. 
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are increased before the first transformation and decreased again after the back 
transformation. 

Figs. 14. 1^ - 14. 241 show this error for the two couphngs rj = 0.2933 and rj = 0.1726. 
At each couphng the solutions computed with a given number of coefficients are 
compared to those, obtained by using twice this number. One can see, that dou- 
bhng the number of coefficients shrinks the error by several orders of magnitude. 
Furthermore at rj = 0.1726 a larger number of coefficients is needed in order to 
keep the error at the same magnitude as for higher couplings. Fig. 14.251 shows 
that keeping the number of Fourier coefficients fixed, the error increases with 
decreasing rj. 

4.4.7 Stability of the DSS Solution 

The stability of DSS solutions might be analyzed in a similar way as the stability 
properties of CSS solutions. Denoting the metric functions and the field by Z{t, z) 
we write the perturbed DSS solution as 

Z(r, z) = ZDss{r, z) + 6Z{t, z), (4.146) 

where Zdss{t, z) is the DSS solution which is periodic in r. The equations are 
then linearized in the perturbations 5Z. The main difference to the corresponding 
problem for CSS solutions is, that the coefficients are not independent of time, 
but depend on time r in a periodic way. Therefore one sets 

5Z{T,z)=e^^5Z{T,z), (4.147) 

where 5Z(r, z) is periodic in r with the period A of the DSS solution. Inserting 
this ansatz into the linearized equations again yields a time-periodic boundary 
value problem (the boundary conditions originating from regularity at the origin 
and the past SSH). This problem can be solved in the same way as above. Due 
to lack of time we had to postpone these computations. 

Instead J. Thornburg adapted the "matrix analysis" fSec. l4.3.'Sj) to perturbations 
of DSS solutions. The method works as described in Sec. 14.3.31 with the only 
difference, that one has to evolve the (slightly perturbed) DSS solution for a 
whole period (or alternatively due to the additional symmetry for half a period) 
instead of integrating only for one time step. The reason for this is, that the 
perturbations depend not only exponentially on time r but also periodically, as 
can be seen from ()4.147|) . In order to extract the eigenvalues A, one has to make 
sure, that the unknown function 5Z drops out of the problem, which is the case, 
when slices are compared, that are half a period apart. 

We mention, that again we expect gauge modes to be detected by this method. 
As for CSS perturbations, there is one gauge mode with A = 1. But now due to 
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Figure 4.19: The error of the subsidiary Einstein equation ()4.145|) for rj = 0.29336. 
The solution was obtained with N = 32. Plotted is the absolute value of the 
expression ()4.145|1 as a function of the spatial variable y at the time steps Xj = 
iA/N for i = 0,4,8, 12, thereby spanning almost half the period. The echoing 
period was computed to be A = 0.5403. 
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Figure 4.20: The same situation as in Fig. 14.191 where this time the solution 
was computed using N = 64. Plotted is the absolute value of the expression 
()4.145|1 as a function of the spatial variable y at the time steps Tj = iA/N for 
i = 0, 8, 16, 24, i.e. at (approximately) the same time steps as in Fig. 14.191 The 
period computed with = 64 is A = 0.5399. From these figures it is clear that 
the error is reduced by several orders of magnitude, when increasing the number 
of Fourier coefficients. 
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Figure 4.21: The error of the subsidiary Einstein equation ()4.145j) for t] = 0.1726. 
The solution was obtained with = 128. Plotted is the absolute value of expres- 
sion ()4.145|) as a function of the spatial variable y at the time steps Tj = iA/N 
for i = 0, 16, 32, 48. The computed value of the period is A = 1.7551. 
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Figure 4.22: The same situation as in Fig. 14. 2H where this time the solution was 
computed using N = 256. Plotted is the absolute value of expression ()4.145|) as a 
function of the spatial variable y at the time steps Tj = iA/N for i = 0, 32, 64, 96. 
With this number of coefficients the period was computed to be A = 1.7521 
Again the error is reduced by several magnitudes, by doubling the number of 
Fourier coefficients. Note that at this small coupling constant [rj = 0.1726) more 
Fourier coefficients are needed to obtain a similar accuracy than in Figs. 14.191 
andlOOlfor r] = 0.2933. 
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Figure 4.23: For rj = 0.29336 the Fourier coefficients of the expression ()4.145|1 
are compared for solutions obtained with = 32 (denoted by "+")and = 64 
(denoted by "x"). Shown is the maximum over the spatial coordinate y = \n{z) 
of the absolute value of the Fourier coefficients plotted against the coefficient 
number. 
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Figure 4.24: For rj = 0.17262 the Fourier coefficients of the expression ()4.145|) are 
compared for solutions obtained with = 128 (denoted by "+") and = 256 
(denoted by "x"). Shown is the maximum over the spatial coordinate y = ln{z) 
of the absolute value of the Fourier coefficients plotted against the coefficient 
number. 
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Figure 4.25: A comparison of the error for a fixed number of coefficients (M = 
32, A'" = 64) for two different values of the couphng constant r] Plotted is the 
maximum over y of the Fourier coefficients of the expression ()4.145|) against the 
coefficient number. Compare these errors to Fig. I4.18[ where the magnitude of 
the coefficients of the variable (p are shown. As the latter increase relative to the 
first coefficients, the solution obtained with a fixed number of coefficients gets 
less accurate. 



the translation invariance in r of the DSS equations, an additional mode^ with 
A = should show up, with the eigenfunction being of the form 

6Z{T,z) = drZDss{r,z). (4.148) 

Again the method suffers from non-convergence with grid resolution. But for a 
certain number of grid points the gauge modes are reproduced rather well. 

Using this method Jonathan Thornburg investigated stability of the DSS solution 
for some values of the coupling constant. He reports [71] that the DSS solution 
has one unstable mode for rj > 0.18. For rj = 0.1726 the numerical results are 
not reliable. 



4.5 The Spectrum of Self-Similar Solutions Rel- 
evant for Type II Critical Collapse 

We summarize the results on existence, properties and stability of self-similar 
solutions in Table I^^Tl As will be explained in Chap. El the stability properties of 

^We note that this mode is not a gauge mode, i.e. it cannot be removed by a coordinate 
transformation (see the footnote in Sec. I4.3.2|) . 
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the solutions studied here will be essential for the dynamics of the model. The 
stable CSS ground state probably is a (singular) end state of time evolution for 
a certain class of initial data. The first CSS excitation and the DSS solution 
having one unstable mode are candidates for intermediate attractors for near 
critical initial data in type II critical collapse. As is clear from Table 14.11 the 
spectrum of these solutions strongly depends on the coupling. For large 77 we 
have constructed the DSS solution, which probably ceases to exist at ~ 0.17. 
The CSS solutions on the other hand only exist for < rj < 0.5. Furthermore 
the existence of marginally trapped surfaces in the analytic continuations of the 
CSS solutions will be relevant (see Chap. Eland the remark in [TT]). 

4.5.1 Comparison of CSS and DSS solutions 

The results of our numerical simulations of type II critical collapse for intermedi- 
ate couplings (described in Sec. 15. 8|) led us to compare the DSS solution with the 
first CSS excitation in the range 0.1726 < 1] < 0.18. In this range of couplings 
both solutions exist. Comparing the profiles of the two solutions one sees, that 
at T] = 0.1726 one can find a DSS phase such that the DSS and CSS solution re- 
semble each other strongly in some fraction of the DSS "backwards light cone"^. 
This resemblance is rather good up to the past SSH of the CSS solution (which 
does not agree with that of the DSS solution). Fig. 14.261 illustrates the situation. 

For rj = 0.1805 on the other hand the agreement (for the "best fitting" DSS 
phase) is not as good as can be seen from Fig. 14.271 

For the following considerations we have to introduce concepts from the theory 
of dynamical systems^, which will also be useful for understanding critical phe- 
nomena (see Sec. 15.111 . Consider the (characteristic) initial value problem for the 
a model in spherical symmetry. As described in Sec. \'2M\ a complete set of initial 
data is given by the field (j) at the initial null surface, 4>o{r) = (f){uo,r). These 
data then are evolved by the means of Eqs. ()2.48p . ()2.55p and ()2.56p . This system 
can be viewed as an infinite dynamical system in the following way: phase space 
is the set of all (asymptotically flat) initial data. An initial configuration (j)o{r) 
thus corresponds to one point in the (infinite dimensional) phase space. Time 
evolution (Eqs. ()2.48p . ()2.55p and ()2.56|) ) of the initial data 0o(^) corresponds to 
a trajectory (an orbit) in phase space. 

In adapted coordinates the CSS solution is independent of time. Time evolution 
maps these data onto themselves, so this solution is a fixed point of the system^. 
An initial configuration that corresponds to the DSS solution, is mapped onto 

^By this sloppy formulation we mean the region on a null slice r = const (u — const) 
bounded by the intersection of the past SSH with this slice. 

^Textbooks for dynamical systems are e.g. ^ and 0]. deals with infinite dimensional 
dynamical systems 

^neglecting the fact, that the CSS solution is not asymptotically flat 
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< < 0.069 


0.069 <7]< 
0.15 


0.15 <7] < 
0.17 


0.17 < 7^ < 0.4 


0.4 < 77 < 0.5 


0.5 < 77 < oo 


CSS 

groundstate 


CSS ground state exists 




regular up to 
future SSH 


analytic extension beyond past SSH contains marginally 
trapped surfaces 




stable 


?? 




CSS 1st 

cXClXdXlOIl 


CSS 1st excitation exists 




regular up to future SSH 


analytic extension beyond past SSH 
contains marginally trapped surfaces 




one unstable mode 


?? 




DSS 
solution 




DSS solution exists 




?? one unstable mode 



Table 4.1: The spectrum of self-similar solutions relevant for type II critical collapse 




Figure 4.26: Comparison of CSS and DSS solution at the coupling rj = 0.1726, the 
lowest at which we explicitly constructed the DSS solution. Plotted are the CSS 
solution (dashed line) and the DSS solution (dots) at a special instant of time, 
where the amplitude of the field (poss is maximal. The vertical lines denote the 
intersections of the slice r = const with the past SSH of the CSS, respectively 
DSS solution. At this coupling both solutions agree rather well - though not 
exactly - up to the past SSH of the CSS solution. 
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Figure 4.27: The same situation as in Fig. 14.261 for a coupling of rj = 0.1805. At 
this coupling one cannot find an instant of time in the DSS solution, for which the 
shape of the field resembles that of the CSS solution as closely as for r] = 0.1726. 
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itself after one period A. The DSS solution therefore can be viewed as a limit 
cycle of the system. 

Our dynamical system depends on a parameter, the coupling constant rj. Ex- 
istence and stability properties of fixed points and limit cycles in a dynamical 
system may depend on such an "external" parameter. In particular the number 
of fixed points and limit cycles might change at some critical value of the parame- 
ter rjc- This "process" is called bifurcation. There are so-called local bifurcations, 
where the appearance of a new fixed point (or limit cycle) is connected to a change 
in stability of the already existing fixed point. And there are global bifurcations, 
where the fixed point keeps its stability properties (for the possibilities of global 
bifurcations in two dimensions see e.g. Chap. 8.4 in |68]). 

Here we are interested in homoclinic loop bifurcations, which are global bifur- 
cations. Fig. 14.281 shows a schematic picture of a homoclinic loop bifurcation 
(for the simple case of phase space being two dimensional): for rj < rjc a fixed 
point with one unstable direction exists. Increasing the parameter towards rjc 
the stable and unstable manifold bend more and more towards each other until 
at rj = rjc they merge and a homoclinic loop develops: one can "leave" the fixed 
point along the unstable manifold and return to it via the stable manifold. Of 
course such a "motion" would take infinite time. For r] > 7]c the homoclinic loop 
separates from the fixed point as a limit cycle. Stable and unstable manifold 
of the fixed point break apart. During this "process" the fixed point does not 
change stability. In principle the emerging limit cycle can be either stable or 
unstable. Approaching the critical value of the parameter from above, it is clear 
that the period of the limit cycle diverges in the limit rj ^ rj'^. For a homoclinic 
loop bifurcation jHE] gives the scaling of the amplitude as 0(1) and of the period 
of the limit cycle as 0(ln(?7 — r^c"))- 

Returning to our situation, we concentrate on the "vanishing" of the DSS solution 
at rj ~ 0.17. We summarize some features of this process: 

1. CSS and DSS solution "come close" in phase space as one approaches rj ~ 
0.17 from above. They lie "farther apart" for bigger rj. 

2. The first CSS excitation does not change stability around rj ~ 0.17. 

3. The DSS period A rises sharply and seems to diverge at ?7 ~ 0.17. 

4. The amplitude of the DSS oscillations is 0(1). 

This suggests, that the DSS solution "emerges" from the CSS solution at r/ ~ 0.17 
in a bifurcation. From 2. one might conclude, that the bifurcation is not a local 
bifurcation (as would be e.g. a Hopf bifurcation) but rather a global one . 3. and 
4. suggest that the bifurcation is a homoclinic loop bifurcation^^ . 

^"^3. and 4. would also fit to an "infinite loop bifurcation", but we consider tliis as unlikely. 
^^The DSS solution has the additional symmetry (/)£)55(t+A/2, z) ~ —4>dss{t, z). Therefore, 
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Figure 4.28: Example of an homoclinic loop bifurcation: ai rj < rjc fixed point 
with one unstable direction exists. At rj = rjc the unstable manifold and the 
stable manifold merge to form a homoclinic loop. For rj > rjc the homoclinic 
loop separates from the fixed point as a limit cycle. The fixed point does not 
change stability throughout. 

Assuming, that we really deal with a homoclinic loop bifurcation at ric — 0.17, 
we can give the following arguments for the behavior of the period A of the DSS 
solution: For rj slightly bigger than rjc, where CSS and DSS are already "close", 
we separate the period into the time T, the DSS solution spends in the vicinity of 
the CSS solution and the remainder Tj-em- As the DSS solution has the additional 
symmetry 0(r + A/2,z) = —(f){T,z), the DSS solution comes close to the CSS 
solution twice (to 4>css and —(pcss) during one period. Therefore we can write 

A^2T + 2Trem- (4.149) 

If DSS is close to CSS we can expand the DSS solution in terms of the CSS 
solution and its perturbations: 

(pDSsir, Z) = (j)CSs{z) + S(f)unstable{r, z) + S(f)stable{r, z). (4.150) 

Note that in this equation the coordinates r and z are adapted to the symmetry of 
the CSS solution, in particular the DSS solution is not periodic in the coordinate 
r. This fact does not matter here, as we are only interested in the local behavior 
in the vicinity of the CSS solution. 

if there is a phase at which the DSS solution resembles the CSS solution (j)css, within the same 
cycle there is another phase (separated by A/2) at which it resembles —(j)css- Strictly speaking 
a heteroclinic loop connecting (pcss to —/pcss forms at the bifurcation point. The bifurcation 
therefore should be called heteroclinic loop bifurcation. Nevertheless we prefer to stick to the 
term homoclinic here, because we think this expresses the essential features. 
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We define ri to be the moment of time, wliere tlie stable modes fiave slirinked 
to order e (||50s4aMe|| = e in some suitable norm). According to tlie definition 
of a fiomoclinic loop bifurcation, tlie admixture of the unstable mode in ()4.150|) 
depends on rj and goes to zero as rj tends to rjc- Therefore we can always find a 
value ?7o such that the norm of the unstable mode at this moment of time ti is less 
than e for all rjc < rj < rjo. We define T2 > ti to be the moment of time, where 
the unstable mode has grown to order e (| |(50unsta6ie| I = e)- From the stability 
analysis we know, that the CSS solution has one unstable mode with eigenvalue 
Ai, which does not depend strongly on 77. Writing ScpunstaUe = ^o^^^^vi^), the 
time T elapsing between ti and T2 is given by 

r= -^lnio + ^In^, (4.151) 
Ai Ai \\y\\ 

where Aq denotes the amphtude of the unstable mode at the time Ti. 

Now the only expression in ()4.15ip that depends on the parameter rj is the am- 
plitude Aq. (We neglect the ?7-dependence of A 1 as A 1 is only slowly varying with 
Tj). By definition it should go to zero for t] rjc- If we assume further that Aq is 
a regular function of rj — rjc, namely Aq^t]) = a{rj — rjc) +0{{rj — rjc)^), we obtain 
the following formula 

T = — — \n{ri — rjc) + const. (4.152) 
Ai 

We may assume further that for 77 close to t]c, the remaining part of the period 
can be approximated by a constant, Trem — const. Therefore we have 

2 

A{r]) = -—\n{r] - r]c) + const. '^^ (4.153) 
Ai 

Fig. 14. Ifil shows the period A fitted against the function f{ri) = —aln{r] — r]c) + b. 
As stated there, the fit gives rjc — 0.17 and a = 0.36278. According to Eq. ()4.153|) 
this would correspond to an unstable eigenvalue Ai = 5.51298. The "true" value 
of Ai at = 0.17, computed with the shooting and matching method as in 
Sec. 14.3. H is Ai = 5.14282. The relative difference of these quantities is ~ 7%. 
This correspondence of numbers gives a strong support to the hypothesis of the 
homoclinic loop bifurcation. 



^^This argument was pointed out by C. Gundlach to us, however on the basis that a 
second unstable mode of the CSS solution appears at the bifurcation. 
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Chapter 5 



Type II Critical Behavior of the 
Self-gravitating a Model 

This chapter finally deals with type II critical phenomena of the self-gravitating 
SU(2) (J model in spherical symmetry. This model has already been studied in 
its limits of strong coupling {jj — > oo) by Liebling [SI] and weak coupling {j] = 0) 
independently by Bizon et al. ITU] and Liebling et al. [S^. They find type II 
critical behavior governed by self-similar solutions in conformity with Table 14.11 
From these results and our knowledge of self-similar solutions (Chap.EJ we expect 
critical phenomena to depend strongly on the coupling. In particular we expect 
the critical solution to change from CSS to DSS in some intermediate regime of 
couplings around ~ 0.17 (see Sec. 15. 8|) . 

As is clear from Chap. |3] (especially Table l^?T|) . for small couplings we have to 
consider the possibility of the formation of naked singularities - according to the 
stable CSS ground state - for a certain class of initial data. This is investigated 
in Sec. 15. HI 

In agreement with Table 14.11 we essentially find three different types of critical 
behavior: for small couplings the critical solution is CSS (see Sec. I5.7|l . while for 
large couplings we have DSS critical behavior (see Sec. I5.(ij) . And for some inter- 
mediate range of couplings 0.15 < ?7 < 0.18 we find that the intermediate asymp- 
totics of near critical evolutions show a behavior which we call "episodic CSS": 
at intermediate times we see a repeated approach to the first CSS excitation. 
These "episodes" are part of an approach to the DSS solution at couplings where 
the latter exists, and still have some resemblance with discrete self-similarity at 
couplings where we think the DSS solution does not exist anymore (see Sec. 15. 8|) . 
With increasing coupling the "CSS episodes" get less pronounced, while at the 
same time the number of episodes (or cycles) increases. To our knowledge this 
sort of transition from CSS to DSS as the critical solution, which is in very good 
agreement with our results obtained by a direct construction of the self-similar 
solutions and the hypothesis of a homoclinic loop bifurcation of Chapter has 
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not been observed in the context of critical phenomena of self-gravitating systems 
up to now. 

The results for large couplings are summarized in . Part of the phenomena in 
the transition region is described in A more complete description is given 
in gn]. 

In order to avoid confusion originating from the inconsistent use of the notion 
"the critical solution" in the literature, we fix our convention as follows: by the 
critical solution we denote the intermediate attractor, which by the means of its 
"stable manifold" separates two different end states in phase space. Solutions that 
approach the critical solution in some intermediate regime of time are called near 
critical solutions (evolutions, data etc.). The member of a family of initial data 
with p = p* is called critical data. With this nomenclature we have constructed 
the respective "critical solutions" in Chapter HI and deal with the evolution of 
"near critical data" in this chapter. It is the aim of the bisection procedure (see 
Sec. 15. 5j) to approximate "critical data" as close as possible, but of course they 
are not realized numerically. 

5.1 Introduction to Critical Phenomena 

There is a couple of excellent review articles on critical phenomena. An elemen- 
tary introduction to critical phenomena is given by Choptuik, who pioneered the 
work on this field, in [20] . Gundlach's reviews (3^1 EZ| give more details as well as 
an overview of the models studied and the phenomena found. We also mention 
the review by Brady and Cai [T^. At the moment the most recent reference lists 
can be found in the article by Wang |22j and on Choptuik's home-page jT7j . 

The field of study can be explained as follows: in their long time evolution iso- 
lated (asymptotically fiat) self-gravitating systems are supposed to evolve to some 
stationary end state, e.g. to a black hole, a stable star or flat space. According 
to this small number of distinct kinds of end states, the space of initial data is 
divided into basins of attraction. The "boundaries" of these basins and their 
"vicinities" are the scope of studies of critical phenomena. 

In the simplest models (e.g. the massless Klein-Gordon field studied by Choptuik 
in spherical symmetry), where only two different end states are possible, namely 
black holes and Minkowski space, "small" initial data, i.e. initial data, that do not 
deviate too much from Minkowski data, will finally disperse to infinity, leaving 
flat space behind, whereas for "strong" initial data, part of the mass present in 
the initial slice will be trapped and a black hole will form. 

Technically one constructs a one parameter family of initial data, parametrized 
by p, such that for large values of p the data evolve to a black hole, whereas 
for small values of p the data disperse. E.g. for the SU(2) a model in spherical 
symmetry, and working on null slices, Eqs. ()2.48|) . ()2.55|) and ()2.56|) show, that 
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Figure 5.1: A family of Gaussian initial data (f)o{r;p) = pr'^exp{—{r — ro)^/o"^), 
where the width a and the center tq are fixed and the amplitude p parameterizes 
the family. Depending on the value of p, the data will finally evolve to a black 
hole or disperse. 

a complete set of initial data is given by the shape of the field (f)o{u = 0,r) 
at the initial null slice. A one parameter family of initial data then can be 
modeled e.g. by a Gaussian with fixed width and center and the amplitude 
serving as the parameter. For a family constructed this way, there will be a value 
of the parameter p, denoted by p*, which separates initial data that disperse 
{sub- critical data with psub < P*) from those, that form a black hole {super- 
critical data with Psuper > P*)- Phenomena that occur for initial data with p c::^ p* 
are called critical phenomena (See Fig. 15. ip . One of the original key questions 
was, whether the black hole mass for slightly super-critical data can be made 
arbitrarily small (such that the black hole mass as a function of the parameter p 
would be continuous) or has a finite value (such that there would be a mass gap). 
The answer is, that depending on the model under investigation, both behaviors 
can occur. In analogy to statistical physics one distinguishes between two types 
of critical phenomena, type I where the black hole mass shows a mass gap, and 
type II, where the black hole mass is a continuous function of p — p*. 

In the following we will concentrate on type II critical phenomena and only refer 
to the other possibility at the end of this section. The first model, for which crit- 
ical phenomena have been investigated, was the self-gravitating massless Klein- 
Gordon field. This was done numerically by M. Choptuik JSJE]. In order to 
resolve all the features, including self-similarity, he had to develop a sophisticated 
numerical algorithm, which refines the numerical grid, when variations occur on 
too small scales to be resolved. Other models have been studied, including e.g. 
gravitational waves in axial symmetry, perfect fluids, the Einstein- Yang-Mills sys- 
tem etc. For the most recent lists of references see j77] and the bibliography on 
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Matt Choptuik's home-page |17j . 

These investigations showed, that the behavior of near critical evolutions can be 
characterized by three main features 

• Scaling, 

• Self-similarity, 

• Universality. 

Scaling relates the black hole mass of super-critical data, as well as other quan- 
tities, that have dimension of length or any power thereof to the parameter p in 
the initial data. One finds, that the black hole mass scales as 

rriBH ~ {p-p*y, (5.1) 

where the exponent 7, called the critical exponent, is independent of the family 
of initial data, although it depends in general on the model. 

The second general feature of type II critical phenomena is that near critical 
evolutions spend their intermediate asymptotics in the vicinity of a self-similar 
solution (which can be either continuous or discrete, depending on the model), 
before they actually decide whether to disperse or to form a black hole. 

The third important point is universality, which means independence of the above 
features of the family of initial data. 

An explanation of these phenomena can be given in the language of dynamical 
systems. Suppose the matter model admits a self-similar solution. For simplicity 
we concentrate on continuously self-similar solutions 4>css{z). Suppose further, 
that the CSS solution has exactly one unstable mode with eigenvalue A and 
that an initial configuration, which corresponds to the CSS solution plus a small 
admixture of the unstable mode leads to either black hole formation or dispersion, 
depending on the overall sign of the perturbation. In the simplest case the "stable 
manifold" of this solution divides the phase space into sub and super-critical data. 
Fig. 15.21 shows a sketch of this scenario using a "phase space picture" . 

Then general near critical initial data are attracted by the CSS solution via the 
stable modes, until they are close to the CSS solution. In this vicinity the solution 
can be written as a small perturbation of the CSS solution 

0(M,r) = (Pcssi-^) + C{p){u* - u)-^y{-^) + 6(^stabUu,r), (5.2) 
u — u u* — u 

where the eigenvalue A is real for all known examples and positive. The amplitude 
of the unstable mode contains information on the initial data, in particular it de- 
pends on the parameter p. For p = p* the unstable mode is tuned out completely. 



105 



black hole 




flat space 



Figure 5.2: A schematic picture of phase space. Every point in this figure cor- 
responds to one configuration 0(r, m = const). In adapted coordinates the CSS 
solution is a fixed point of the system, it therefore is drawn as a point (large solid 
circle). The CSS solution has one unstable mode, the "stable manifold" therefore 
is of co-dimension one. The straight line at the left of the figure represents a one 
parameter family of initial data, with parameter p. The value p* corresponds to 
those initial data, that "start out" on the "stable manifold" and are completely 
attracted to the CSS solution (where they arrive only in the hmit r — > oo) . For 
p > p* the configuration is initially attracted by the CSS solution via the stable 
modes until the (initially very small) admixture of the unstable mode takes over 
and pushes the solution towards black hole formation. For p < p* the final state 
is fiat space. Any one parameter family of initial data, cutting the "stable man- 
ifold", will show the same near critical phenomena. Two remarks are in order 
here: a DSS solution, being periodic in the adapted time r, corresponds to a limit 
cycle, and should therefore be drawn as a cycle, with near critical data spiraling 
in and out. Second, this sketch does not claim, that the "stable manifold" is in- 
deed a manifold. It is just a very helpful abstract picture, modeling the essential 
facts, that have been observed. 
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and the configuration evolves towards the CSS solution, therefore C{p*) = 0. For 
near critical data, we have C{p) = ^{p*){p — p*) + {0{{p — p*Y))- 

From ()5.2|) we can estimate the time for which the solution stays in the vicinity of 
the CSS solution. Fix ti to be some instant of time, where near critical solutions 
are already in the vicinity of the CSS solution, and therefore fl5.2|) is valid. Let 
furthermore T2 = ti+T he the instant of time when the amplitude of the unstable 
mode has grown to be e, then the time T spent in the vicinity of the CSS solution 
is given by 

6 = Cip)e'^, (5.3) 

with C{p) = C{p)e^^\ or 

^ 1 , const . , 
T = -ln . 5.4 

So a near critical solution spends longer and longer time in the vicinity of the CSS 
solution, when p comes closer and closer to p*, until for p = p*, the logarithmic 
time goes to infinity. 

In order to explain the scaling of the black hole mass, we "redefine" our family 
of initial data in the following way: we fix the time Up, where the stable modes 
are already negligible compared to the unstable mode, and the amplitude of the 
unstable mode has grown to be e 

e = C{p){u*-Upr\ (5.5) 
e = ^{p*){p-p*){u*-up)-' (5.6) 



So 



dp 



or resolved for u* — Up 



u* — Up = const{p — p*) ' (5.7) 

where the constant contains e and some information on the original family of 
initial data, but is independent of p. This way we have constructed a family of 
initial daXa. (pQir^p) = (pcss{f^ / {u* —Up))+ey{r / {u* —Up)) , which depends on r only 
via the ratio r/{u* — Up). As the field equations are scale invariant, a solution 
(j){u^r) with initial conditions 0o(^) implies the existence of a one parameter 
family of solutions 0cr(M, r) = (j){au, err) with initial conditions {(pa)o{^) = (po{o'r). 
From this it follows, that the evolution of our one-parameter family of initial data 
gives a one-parameter family of solutions of the form 

(/),(«, r) = , ). (5.8) 

U —Up U — Up 

This is valid for the whole future evolution of the data, even when the linearity 
assumptions break down. 
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Assume now, that for super-critical data the solution (f){u, r) has an apparent 
horizon at thIu). If we fix z = r/{u* — u), then the apparent horizon is located 
at {uH,rH = z {u* — uh)) and so the metric function j3{r/{u* — u),u) diverges 
when u — > uh- Therefore the rescaled solutions (3p{r/{u* — u),u) = I3{r/{u* — 
u),u/ {u* — Up)) diverge if m ^ {u* — Up)uH, or r ^ {u* — Up)rH- The mass of the 
apparent horizon m{z;p) = rH{z;p)/2, measured at constant z, therefore scales 
as {u* — Up) . 

Remark: the above analysis, of writing near critical data as the CSS solution 
plus a perturbation is valid only up to some finite radius. Outside this region the 
near critical solutions - being asymptotically flat - will differ considerably from 
a small perturbation of the CSS solution. Nevertheless, if the region outside does 
not influence the black hole mass the latter scales as 

"m-BH ~ {u* — Up) = const{p — . (5.9) 

So we have derived the scaling law ()5.1|) with the additional information, that 
the critical exponent 7 is related to the eigenvalue of the unstable mode via 

7 = ^. (5.10) 

This relation was derived independently by Koike et al. jlZj and Maison 
(Evans and Coleman [23] first suggested to look at the linear stability of the CSS 
solution in order to get an estimate on the critical exponent 7.) Another well 
defined quantity as pointed out by Garfinkle and Duncan 27 , which exhibits 
scaling, is the maximum of the Ricci scalar at the axis for sub-critical data, 
the maximum taken over a whole evolution, max7^(M,0). As this quantity has 

u 

dimension of 1/ length"^, it should scale as 

max7^(M, 0) ~ {p* - p)'^/^. (5.11) 



If a DSS solution is the critical solution, the scaling law undergoes some modifi- 
cation, in that a small wiggle is overlaid. The derivation, as first given indepen- 
dently by Gundlach [32] and Hod and Piran [Hj, is analogous to the CSS case, 
a first difference arising in ()5.2|1 . where now the DSS solution and its unstable 
mode have an additional periodic dependence on r. The family of initial data 
constructed as above, therefore depends periodically on the parameter r^. 

Again the equations are scale invariant, therefore the solutions to initial condi- 
tions 0o('";p) behave as 
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Any quantity of dimension length should therefore scale as e~'^p({Tp)f{Tp) = 
e~'^pf{Tp), where /, C and / are periodic functions of their argument Tp. So we 
have 

mBH = ci{p-p*y/^f{-\nci - -ln(p-p*)) (5.14) 

or ^ ^ 

InniBH = Inci + - \n{p - p*) + /(- Inci - - \n{p - p*)), (5.15) 
A A 

where / = ln(/). Now / is periodic in ln{p — p*) with period A/7, or since 
the metric functions have the additional symmetry of consisting only of even 
frequencies, and therefore have a period of A/2, the period of / is rather half 
this value, i.e. Note that there is only one constant ci, which depends on the 

family. The scaling exponent 7 = 1/A and the periodic function / are universal. 

We close this section with some words on type I critical phenomena, which have 
been observed in several models, e.g. in the Einstein- Yang-Mills system j2Tj or 
the massive minimally coupled scalar field In this type of transition the 

intermediate asymptotics is governed by an unstable (the stable manifold being 
of co-dimension one) static solution or a solution that is oscillating in time. Again 
the "life-time of the critical solution", that is the time a near critical solution 
spends in the vicinity of the intermediate attractor, scales according to ()5.4|1 . 
The major difference to a type II collapse is that the black hole mass for slightly 
super-critical data is finite, i.e. the black hole mass as a function of the parameter 
p — p* is discontinuous. The magnitude of the mass gap or the fraction of the 
mass of the intermediate attractor that is radiated away by shghtly super-critical 
data depends on the model. 



5.2 Limits of Weak and Strong Coupling 

The self-gravitating SU(2) a model has already been investigated with respect to 
critical phenomena for the two extremes of coupling, rj = and 77 — 00. 

The case t] = corresponds to the SU(2) a model on fixed Minkowski background 
and was investigated independently by Bizon et al. jlO and Liebling et al. 
They looked at the threshold between dispersion and blow up (of the first deriva- 
tive of with respect to r) at the origin, which is governed by the CSS ground 
state. They found that the solution at the threshold is the first CSS excitation. 
A quantity that can be used to examine the scaling, is the maximum (over time) 
of the energy density at the origin. It was found, that this quantity shows scaling 
with an exponent 7 = l/Xcss, where Xcss is the eigenvalue of the unstable mode 
of the CSS solution for 17 = ffl. 

The case — > 00 on the other hand corresponds to the self-gravitating a model 
with three dimensional flat target manifold (]R^,(5ab). This is natural, as the 
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coupling r] corresponds to the inverse of the curvature of the target manifold, the 
curvature goes to zero, as 1] tends to infinity. It is also easy to see, by e.g. looking 
at Eqs. ()2.48|) . ()2.55|) and ()2.56|) . Defining (p = ^ 0, rewriting the equations with 
respect to and taking the limit r] oo gives the following system of equations 

r 

□0 = ^, (5.16) 

which are precisely the equations for the a model with fiat target manifold in the 
hedgehog ansatz. 

This model was investigated by Liebling [ST], where in addition he considered 
a potential. As was explained in Sec. I4.1.3| such a potential is asymptotically 
irrelevant, the critical behavior of the models with and without this potential 
should therefore be the same. Liebling found that the critical solution at the 
threshold of black hole formation is DSS with an echoing period A = 0.46 and 
a scaling exponent 7 = 0.119. Note, that the value of the echoing period nicely 
fits to the period of the DSS solution aX rj = 100, described in Sec. 14.4.51 

From these rather different critical phenomena at the limits of very small and 
very large couplings, one can infer that there will be a transition of the critical 
solution from CSS to DSS as the coupling is increased. 



5.3 Possible End States 

As is clear from the last section, criticality is only defined with respect to the 
end states. Usually these two different end states would be black hole formation 
and dispersion to infinity. If the model allows also for other stable stationary 
configurations, e.g. for stable static solutions, then these can be considered as 
possible end states as well. 

As was already mentioned in Sec. 12.1. If the SU(2) a model does not allow for 
static asymptotically flat solutions. Therefore it is natural to investigate the 
transition between black hole formation and dispersion. 

On the other hand, as we have seen in Chap. 01 the spectrum of CSS solutions 
contains a stable ground state. Furthermore on fixed Minkowski background this 
stable ground state governs the "long time behavior" of strong initial data, as was 
described in jHH and the intermediate asymptotics of near critical data between 
this singularity formation and dispersion is ruled by the first excitation of this 
CSS family. 
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As the ground state and its stability properties persist, when gravity is turned on, 
it is reasonable to expect, that a naked singularity is a possible "end state" for 
strong (but not too strong) initial data, at least for couplings less than rj ~ 0.069. 

As a first step we tested, whether "semi-strong" (but otherwise arbitrary) data 
really develop towards the ground state. For r] = 0.01 and t] = 0.05 we evolved a 
Gaussian ()5.17|1 with width a = 1.0 and center ro = 5.0 and chose an amplitude 
p = 0.03, which is neither too big (such that there is no black hole formation) 
nor too small (such that the data don't disperse), but otherwise arbitrary. What 
we find is, that for such data (chosen rather arbitrarily) the solution evolves 
towards the CSS ground state (see Figs. 15.31 and 15. 5|) and stays there, until the 
lack of resolution near the origin causes the numerical data to break away. (The 
numerical evolution then does not represent a solution to the Einstein equations 
anymore). Furthermore we find, that — stays away from 1 during the evolution, 
(see Figs 15.41 and 15. 6|) We can conclude from this, that indeed a naked CSS- 
singularity is the generic end state of "intermediately strong" initial data. 

Due to numerical difficulties, we were not able to determine the value of the 
coupling constant, from which on black hole formation is the only possible end 
state for "strong" initial data. Presumably it lies close to 1]^ ~ 0.069, where 
the CSS ground state has marginally trapped surfaces. We can only state, that 
for 1] > 0.09 we detect black holes as super-critical end states, which show the 
expected scaling, and that for 77 > 0.1 the black hole masses show second order 
convergence. (For details see Sec. 15. 7|) . 

5.4 The first CSS excitation for 0.15 < 77 < 0.5 

In Section 14.2.31 we demonstrated, that the first CSS excitation, if continued 
analytically beyond the past SSH, contains marginally trapped surfaces for t] > 
0.152. This fact might prevent this solution to play the role of a critical solution 
between dispersion and black hole formation for couplings 0.15 ^ rj < 0.5. 

In order to investigate this, we matched a certain class of asymptotically flat 
data^ to the first CSS solution, the matching point being the past SSH and the 
matching condition being such that the resulting data were C^. J. Thornburg 
evolved these data (with several values for the parameters) numerically for t] = 0.2 
and found, that they developed an apparent horizon outside the past SSH. It is 
reasonable to assume, that these data will show the same behavior for couplings 
0.2 < f] < 0.5. For couplings 0.15 < rj < 0.2 the numerical evolution does not 

^ These data are given in the following way: at the past SSH the CSS solution is matched to 
a cubic polynomial such that the resulting data are . At some distance away from the past 
SSH the cubic polynomial is matched to a Gaussian, the matching again being C^. The two free 
parameters for these data are the location of the second matching point and the width of the 
Gaussian, all the other parameters of these data are used to achieve the required smoothness 
at the matching points. 
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Figure 5.3: Late time behavior of initial data as described in the text for rj = 0.1. 
This figure shows the evolved field (lines-points) - moving from right to left as 
time proceeds - for several time steps. The solution clearly comes close to the 
CSS ground state (lines), and stays there until the grid resolution at the origin 
becomes too sparse and the evolved data break away from a solution of Einstein 
equations (this last time step shown in the plot is already "after" the culmination 
time u*). The culmination time u* of the CSS solution was determined by the fit 
of a single time step to be u* = 11.346. For the other time steps the CSS solution 
was shifted according to r = z{u* — u). The past SSH of the CSS ground state 
is located where the field equals tt/2 ~ 1.58. This shows that the region, where 
the evolved data and the CSS ground state agree, extends some way outside the 
past SSH. 
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Figure 5.4: The same scenario as in Fig. 15. 3| where this time — is plotted. Note 
that — is far from being unity everywhere in the evolved grid. This shows, that 
the singularity, which is approached via the CSS ground state, is in general not 
shielded by an apparent horizon. 
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Figure 5.5: For rj = 0.05 the same initial data as in Figs. 15.31 and 15.41 are evolved. 
This figure shows the evolved field (lines-points) for several time steps between 
u = 10.1448 and u = 11.2416. Again the solution clearly comes close to the 
CSS ground state (lines) and breaks away due to insufficient resolution near the 
origin, (clearly the "latest" time step plotted suffers from insufficient resolution) 
The culmination time u* is determined via the fit to be u* = 11.24085. 



0.45 




r 



Figure 5.6: Same data as in Fig. 15.51 This time ^ is plotted. Note that 
although growing in time, ^ does not come close to one anywhere in the slice 
before the culmination time is reached. This means that these initial data lead 
to the formation of a naked singularity. 
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yield definite results, because at these couplings it is easy to construct data, which 
do not form an apparent horizon before the evolved data break away from the 
CSS solution due to numerical errors, which correspond to an excitation of the 
unstable mode. If one could eliminate (or diminish) these numerical errors, it is 
likely, that alos for these couplings all data of the above described class would 
form an apparent horizon. At this point we cannot decide this. 

If the behavior observed for rj = 0.2 is generic, it is clear that the CSS solution 
cannot be found by the means of a critical search between dispersion and black 
hole formation, because the solution itself (matched to asymptotically flat data) 
evolves to a black hole - and so would small perturbations independently of the 
sign of the admixture of the unstable mode. In other words, the CSS solution 
does not lie at the boundary between dispersion and black hole formation. 

Indeed, as stated in Sec. 15.61 the first CSS solution does not show up as a critical 
solution for t] > 0.2. In the transition region 0.15 ^ f] ^ 0.18 as described in 
Sec. 15. 8| the CSS solution appears in the "CSS episodes" of near critical evolu- 
tions, but clearly it is not "the critical solution", i.e. the intermediate attractor, 
whose unstable mode is tuned out by bisection. 



5.5 Critical Searches - Setup and Extraction of 
Results 

In order to investigate critical behavior we used several families of initial data, 
namely a "Gaussian" 

0o(r) =p rV^'^-'^'')'/"', (5.17) 

with the center ro usually set to be 5.0, the width a fixed to be either 1.0 or 2.0 
and the amplitude p serving as the parameter. For large couplings also the family 

Mr) = -4Ar2 (^^^^^ e-("-''«)'/p', (5.18) 

generated by keeping the center rg = 5.0 and the amplitude A fixed and varying 
the width p as the parameter was used (see ^^)■ 

Furthermore we tried a "double Gaussian", 

0o(r) = p ^2g-(r-ro)V-^ + A (.-^-^-r^^^ I -I ^ (5.19) 

with the second Gaussian fixed (A = 0.001, r2 = 7.0, = 0.5), width and center 
of the first Gaussian fixed, a = 1.0, tq = 5.0 and p being the parameter. 

For a fixed value of the parameter the initial data 0o were evolved using the DICE 
code (see App. until for large couplings either a black hole formed (for the 
numerical criterion for a black hole see App. [n]) or the field dispersed (most of 
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the times measured via max2m/r less than some small value, e.g. 10 ^, which 

r 

might of course depend on the coupling; see also remark below). For very small 
couphngs, where we expected naked singularities as super-critical end states, an 
evolution was defined to be super-critical, whenever the errors grew above some 
limit ^. 

Starting with some value p G [pmin,Pmax], the interval chosen such that pmin leads 
to dispersion, while Pmax leads to a super-critical end state, the parameter p was 
driven towards p* by bisection: for pi > p*, the interval [pmin,Pi] was halved to 
give the new value of the parameter p2 = {pmin +Pi)/2, at the same time the 
interval was reset to \pmin,Pmax = Pi] and so on. (This description applies if 
increasing the parameter makes the initial data stronger, as is the case e.g. for 
p being the amplitude of a Gaussian. For p being e.g. the width of a Gaussian 
the parameter has to be decreased to make the initial data stronger.) Such a 
bisection search is limited by floating point errors, so a critical search flnished, 
when {pn - Pn-i)/pn-i < 10"^^. 

Given the result of a critical search, the critical value of the parameter p* was 
approximated by p* = {psub + Psuper) 1'^^ where psuh was the biggest sub-critical 
and Psuper the smallest super-critical value obtained. 

For large couplings, where the critical solution is DSS, the echoing period A and 
the culmination time u* were determined simultaneously from (max— )(n). This 

function of time reflects the periodicity of the DSS solution in logarithmic time 
T = — ln(n* — u) (see Fig. 15. 7j) . A perl script, written by Jonathan Thornburg, 
extracted A/2 and u* using the minima of (max— )(u), at times Ui. The times 

5un elapsing between the two adjacent minima at u„ and Un+i are given by 5un = 
g-(n-i)A/2^^^^ a least squares flt of ln(5M„) to the straight line — (n — 1) A -|- const 
gives the echoing period A/2. Furthermore for an exact DSS solution these times 
5un sum up in a geometric series to give u* = Un + 5m„/(1 — e^'^) (for any n), 
from which u* can be calculated. 

In order to examine the scaling of the black hole mass for super-critical data and 
of the Ricci scalar at the axis for sub-critical data, a whole series of time evo- 
lutions was done, starting close to p* and increasing (decreasing) the parameter 
to log(|p — p*\) ~ —10, with steps equally spaced in ln(|p — p*\)- For the details 
of measuring the black hole mass see App. O To extract the scaling exponent 7 
from the black hole masses, a perl script (written by Jonathan Thornburg) least 
squares fltted \n.mBH{x) to the straight line '-yx + k with x = ln{p — p*). For an 
extraction of 7 from the scaling of the Ricci scalar usually In mBH{x) was fltted 
to the straight line by naked eye. 

^Note: it is always a great pleasure to declare the limitations in accuracy of a numerical 
code a "physical state" . 
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5.6 Critical Phenomena for Large Couplings 



For large couplings 0.2 < rj < 100, we find that the critical solution at the 
threshold of black hole formation is discretely self-similar (See t45|). Figs. 15.71 - 
15. Ill illustrate this for = 100. All the runs for these figures were done with the 
family ()5.17|) (width a = 1.0), with router = 30.0 as the initial spatial extension 
of the grid and with = 2000 grid points initially. 

Fig. 15.71 shows max2m/r for rj = 100, which if plotted vs. — ln(u* — u) = r 

r 

is a periodic function with period A/2. (As the metric functions j3 and ^ are 
periodic with period A/2, max2m/r shows the same periodicity.) Fig. 15. 81 shows 

r 

the Ricci scalar at the origin r = 0, Tl{u,r = 0), which behaves as e^'^i?(r), 
where R denotes a periodic function of its argument. Figs. 15.9) and 15.101 show the 
scaling of the black hole mass for super-critical initial data and max7^(u, r = 0) 

u 

for sub-critical data respectively. Fig. 15.91 also shows the superimposed wiggles in 
the mass scaling. Finally Figs. 15.111 and 15.121 show that a near critical evolution 
comes close to the DSS solution at intermediate times for rj = 100 and t] = 0.2933 
respectively. Compared are the field as evolved from near critical initial data, 
and the DSS solution itself, constructed as described in Sec. 14.41 

Table IHTD (taken from gives the scaling exponent 7, the echoing period A/2 
and the value of the critical parameter for families ()5.17j) and ()5.18|1 for various 
values of the coupling constant rj. As can be seen from this table, for very large 
couplings rj = 100, the scaling exponent 7 ~ 0.1185 and the echoing period 
A ~ 0.461 are in very good agreement with the results for t] = 00 reported in 
[HI] (7 = 0.119, A = 0.46). When decreasing the coupling, the scaling exponent 
7 hardly changes (the variation is at most 5%) whereas the echoing period A 
starts to increase at lower couplings, which is in good agreement with the results 
presented in Sec. 14.4.51 Indeed, as can be seen from Fig. 15.131 the echoing period 
Acrit of near critical evolutions matches the period A^ss of the "exact" DSS 
solution, as described in Sec. 14.4.,'Tl for rj > 0.2. 

For T] = 0.18, although the DSS solution still exists, the observed DSS period- 
icity for near critical evolutions is only approximate. This phenomenon will be 
described in more detail in Sec. 15.81 

5.7 Critical Phenomena for Small Couplings 

For small couplings we investigated critical phenomena between dispersion and 
singularity formation {0 < rj < 0.07) and between dispersion and black hole 
formation (0.07 ^ rj < 0.14). We find that the critical phenomena in this range 
of couplings is governed by the first CSS excitation described in Sec. 14.21 
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— ln{u* — u) 



Figure 5.7: max — as a function of r = \n.{u* — u). The echoing period A was 

computed to be 0.4599. This is the easiest way to extract the echoing period 
from a near critical solution. 




Figure 5.8: The Ricci scalar at the center of spherical symmetry as a function 
of r = — ln(M* — u). As was discussed in Sec. 14.1.21 the Ricci scalar behaves like 
7^(r, z) = e~^^i?(r, z), with R being periodic in r. (Remark: the vertical lines in 
the middle of the figure and near the right end are errors, that occur at each grid 
refinement, but which do not seem to influence the time evolution.) 
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Figure 5.9: The scaling of the black hole mass niBH- Plotted are the masses 
obtained from a series of time evolutions (dots) together with the straight line 
/(ln(p — p*)) = j\n{p — p*) + k, where 7 = 0.1189, versus the left axis. At 
the right axis the difference of these functions is plotted. This reveals the "fine 
structure" of the mass scaling, oscillations with period A/27. 
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Figure 5.10: The scaling of max7^(M,0) for sub critical evolutions as a function 

u 

of ln(]9* —p). The overlaid straight line has slope —27. (Remark: the "escaping" 
points in the middle of the graph presumably stem from a grid refinement.) 
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Figure 5.11: Snapshots of the critical solution (pcrit (solid lines) ai t] = 100 compared to the DSS solution ipDSS (dots). The 
snapshots are taken at times = iA/N for N = 16 and i = 0, 2, 4, 6, 8, 10, 12, 14, i.e. spanning one period. 




Figure 5.12: Snapshots of the critical solution (pcrit (solid lines) at = 0.2933 compared to the DSS solution (poss (dots). 
Snapshots are taken at times Tj = iA/N for N — 64: and i = 0, 16, 32, 48, 64. Note that for r = A the critical solution retains 
its shape at r = 0, but is shifted in ln(r) by —A = —0.5399. 
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Figure 5.13: The echoing period Acrit (large sohd dots) as a function of the couphng constant rj, compared to the period 
Adss of the DSS solution (see Fig. I4.15|) . (The values of Acrit are taken from [13].) 
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Table 5.1: This table shows two families of near-critical initial data parameters 
for yarious coupling constants 77. For the Gaussian-like initial data family 15.171 
we use the 'amplitude" A as the parameter p (at a fixed 'width" a = 1), with a 
numerical grid of 16 000 grid points. For the family 15". 18| we use the 'width" a as 
the parameter p (with different 'amplitudes" A for different coupling constants), 
with 8000 grid points. For each coupling constant and each family, the table also 
shows the max 2m jr echoing period A/2 of the near-critical eyolution, and the 
mass-scaling-law critical exponent 7 determined for the entire critical search. For 
T] = 0.18 the DSS symmetry is only approximate (see Sec. 15. 81 for details). This 
table is taken from jlH]. AH the runs for this table were done by J. Thornburg. 

For the couplings where black holes form (0.07 < rj) the black hole mass scales 
according to ()5.H) with a scaling exponent 7 that corresponds to 1/Xcss, the 
relatiye error being at most 3 %. 

For couplings where the end state is the CSS ground state (77 < 0.07) the quantity 
max7^(M,0) for sub-critical data exhibits scaling according to (jSHI}, although 

u 

only for small yalues of {p — p*)- Furthermore this time the scaling exponent 7 
differs from the theoretical prediction 1/Xcss by up to 15 %. The reason for this 
inaccuracy is not clear to us at the moment. One possible reason could be, that for 
these small couplings ln{{p — p*)) / p* < 10~^^, which limits the determination of 
p* due to floating point round off errors, does not allow to reduce the admixture of 
the unstable mode of the CSS solution in the initial data as much as is the case for 
large couplings, where we obserye a beautiful scahng law fSec. 15. 6|) . Whether this 
is indeed the reason could be checked by switching to higher numerical precision. 
Another possible reason might be, that the code does not work as accurately for 
small couplings as it does for large couplings. Conyergence tests would be a first 
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check of this. 



Figs. 15. 1^ - 15.201 ilhistrate the critical phenomena for the couphngs rj = 0.11 and 
7] = 0.1. Figs. 15.141 and 15.161 show the scahng of the black hole mass for super- 
critical data and of the quantity max7^(M, 0) for sub-critical data. Both figures 

u 

show the results for various initial grid resolutions (A^ = 500,1000,2000,4000). 
As described in j3H] the critical value of the parameter p* depends on the grid 
resolution. This was taken into account in Figs. 15. 111 - 15. 161 i.e. for each resolution 
N the corresponding value p*{N) was used. In |45j it was shown, that for large 
values of the coupling the critical value p*{N) shows second order convergence 
with the grid resolution. We note that this is also the case for small (presumably 
V ~ 0.07) couplings. For r] = 0.11 and the initial data family ()5.17|1 with fixed 
width (7 = 1 and variable amplitude the respective values are 

p*(500) = 0.0214965187393766, p*(1000) = 0.0214974987387296, 
/(2000) = 0.0214977388532529, p*(4000) = 0.0214977981419972. 

The differences are 

9.79999352997835 10"°^ 
2.40114523302609 10"°^ 
5.92887442994738 10'^^, 



5pi 


= p*(1000) 


-p*(5000) 


SP2 


= p*(2000) 


-p*(1000) 


SP3 


= /(4000) 


-p*(2000) 



(5.20) 



and the ratios thereof 



Spi 
Sp2 



4.08138308136726, 



^ = 4.0499175035613. 

dps 



(5.21) 



In addition to the convergence of p* (N) we demonstrate convergence of the black 
hole mass. Fig. 15.151 shows the differences of the black hole masses mBH{N) 
of Fig. 15.141 The difference between runs with and 2N grid points, and the 
difference between runs with 2N and 4N grid points multiplied by a factor of 
four are close throughout the shown interval of p — p*. 

Finally Figs. 15.171 - 15.201 deal with the intermediate asymptotics of near critical 
data for t] = 0.1. Fig. 15.171 shows, that near critical data evolve towards the 
first CSS excitation, stay there for some time and then deviate again. Fig. 15.181 
investigates the deviation of the evolved field from the CSS solution. Accord- 
ing to the theory this deviation should be dominated by the unstable mode of 
the CSS solution. An additional complication arises from switching to adapted 
coordinates, in that an error in determining the culmination time u* brings the 
gauge mode (Sec. 14. 3^^ into the game. In Fig. I5.18l we have taken this fact into 
account in the following way: according to the theory the evolved field should 
behave as 0(r, z) = (pcss^z) + e'^'^ y{z) , where we assume that all the stable modes 
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Figure 5.14: The mass scaling for t] = 0.11 for the family ()5.17|) (with fixed 
width (7 = 1 and variable amplitude). Different symbols denote different grid 
resolutions: denote = 500 grid points initially, "x" = 1000, 

= 2000 and "□" N = 4000. The straight line has a slope of 7 ~ 0.185, which 
is in good agreement with 1/Ai = 0.181. (The relative error is 2%.) 
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Figure 5.15: The differences of black hole masses of Fig. l5.l41 for different initial 
grid resolutions: denotes the difference (ms_H")iooo ~ {^bh)2ooo, "x" denotes 
four times the difference (mB//)2ooo ~ ("^5/^)4000- These quantities almost lie on 
top of each other, which shows second order convergence of the black hole mass 
with the grid resolution (see Sec. l(I5j) . 
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Figure 5.16: The quantity ln(max7?.(ti, 0)) for sub-critical data of the same family 

u 

as in Fig. l5.14l and different grid resolutions: "+" denote = 500, "x" = 1000, 
= 2000 and "□" N = 4000. The straight line has a slope of -27 with 
7 ~ 0.185. 




Figure 5.17: The field of a near critical evolution (dots) at intermediate times 
between u = 12.723 and u = 15.150 (moving from right to left). The field 
approaches the CSS solution (solid lines). The culmination time u* (determined 
from the last but second snapshot) is u* = 15.167. The past SSH of the CSS 
solution is located where the CSS solution attains the value 7i/2 for the second 
time. 
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Figure 5.18: The deviation of the critical solution from the CSS solution for 
several time steps, r = 3.4169, 3.539, 3.6069, 3.6791, 3.756, 3.838. Plotted is 
the difference between the evolved and the CSS solution 4>css, i-e. 50(t, z) = 
4>{t,z) - 4>css{z) (dots). Overlaid is the function a(r)(?/„„st(^) + b{T)ygauge{z)) 
(solid lines), where the parameters a(r) and 6(r) were fitted with bare eye such 
that the maxima of the two functions agreed. The gauge mode has to be taken 
into account, because of the error in determining u*. These parameters should 
depend on time r as a(r) = e'^'^ and 6(r) = bo e^^^^'^'^. So from the fit the 
eigenvalue A could be determined. From Fig. 15.191 one sees, that the value of A 
computed from a(r) is close to the "exact" value of 5.5846. Unfortunately 6(r) 
is not so well behaved (see Fig. I5.20|) . 
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Figure 5.19: The parameter a of Fig. 15.191 in dependence of r (dots). The sohd 
hne has a slope of 5.5846, corresponding to the eigenvalue A computed from the 
perturbation analysis. 
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Figure 5.20: The parameter h of Fig. 15.181 in dependence of r (dots). The solid 
line has a slope of —3.3, which would correspond to an eigenvalue A = 4.3. 
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already have damped out. In adapted coordinates with shghtly different u* this 
would read 4>{t, z) = (pcssi^) + ^^^y{.z) + e'^ygauge{z). Therefore we try to fit 
S(j){T,z) = (I){t,z) - (pcssiz) to the function /(r,z) = a{T){y{z) + h{T)ygauge{z)) 
adjusting a and h such that the maxima agree. The fitted parameters then should 
behave as a(r) = a^e^'^ and 6(r) = h^e'^^"^'''^ . Figs. EUHl and 1^201 show that a(r) 
is in good agreement with the above formulae, whereas 6(r) is slightly off. 

5.8 Critical Phenomena for Intermediate Cou- 
plings — Transition from CSS to DSS 

Finally we describe the region of couplings 0.14 < 0.18, where the transition 
from CSS to DSS as critical solution takes place. We know from the last sections 
(Sec. 15. (il and Sec. 15. 7j) that for r/ < 0.14 the critical solution is the first CSS exci- 
tation whereas for rj > 0.2 the critical solution is DSS. Furthermore in Sec. 14.5.11 
we proposed the hypothesis, that the DSS solution merges with the CSS solution 
at ?7 ~ 0.17 in a homoclinic loop bifurcation and does not exist for smaller r]. 
While our numerical results for the stability of the DSS solution at rj = 0.1726 
are not conclusive (Sec. 14.4. 7|1 it is reasonable to assume that the DSS solution 
does not change stability. 

We start by describing the intermediate asymptotics of near critical evolutions. 
In the whole range of couplings we find a behavior, which we call "episodic CSS", 
that is: the near critical solution approaches the CSS solution (j^css: goes away, 
approaches its negative —(pcss, goes away etc. until after a small number of such 
episodes it finally parts to form either a black hole or to disperse. The culmination 
times u* associated which each episode increase with the episodes. 

In the following we will describe the critical phenomena we find for the two 
coupling constants r] = 0.1726, where we have constructed the DSS solution, and 
7] = 0.16, where we think, that the DSS solution does not exist. 

For T] = 0.1726 the evolution could be compared to the DSS solution. We find that 
the DSS solution is approached better and better during the time the solution 
stays in the neighborhood of the "critical hyper-surface", although not as good 
as at higher couplings e.g. at rj = 0.2. In pEl we will quantitatively give the 
"distance" of the near critical solution to the DSS solution in some norm^ for 
various coupling constants. These investigations show clearly that for rj = 0.2 
the near critical data approach the DSS solution quickly (roughly within one 
cycle), stay in the vicinity (with a distance in the above norm of ~ 10^^) before 
they deviate. For rj = 0.1726 on the other hand the approach to the DSS solution 

■^The discretized version of / dr\(j>{uo,r) — <l>DSs{uo,r)\'^ / J dr, where the DSS 



solution is taken at a time, that minimizes this error and is appropriately shifted in Inr and 

' max 

is mm{router(uo),rssH{u)). 
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absorbs all the time the solution stays in the neighborhood and the field only 
comes as close as ~ 10~^, before it deviates. 

A possible reason for this slow approach (as compared to larger couplings) could 
be that the stable modes of the DSS solution damp out much more slowly at 
T] = 0.1726 than at 77 = 0.2. It would be very interesting to test this behavior 
further, in using a higher numerical precision (quartic precision). This way the 
admixture of the unstable mode in the initial data could be reduced, which would 
prolong the time the solution stays close to the DSS solution. It should then be 
possible to observe a closer approach to the DSS solution at rj = 0.1726. 

Combining the approach to the DSS solution with the fact that CSS and DSS lie 
"close" , it is clear that we observe the above described CSS episodes. 

For T] = 0.16 the CSS episodes are illustrated in Fig. 15.211 We also studied 
the way the evolved solution deviates from the CSS solution at the last but 
one episode. Fig. 15.221 shows the deviation 50 together with the fitted functions 
f{T,z) = a{T){yunstabie{z) + b{T)ygauge{z)) defin ed as in Sec. EHI Clearly th e fits 
are not as good as for 77 = 0.1 (see Fig. I5.18|) . Nevertheless the fits (Fig. 15.231 
and I5.24|) show that the maximum grows exponentially according to the unstable 
mode of the CSS solution. 

The explanation for the episodes in this case is not as straightforward as for 
7] = 0.1726. A possible explanation would be, that although the DSS solution 
does not exist, there are still orbits in phase space, which "mimic" a DSS critical 
solution, i.e. orbits, that do not close exactly, but nevertheless act as intermediate 
attractors. 

Taking this, one would expect, that the black hole mass as well as the scalar 
curvature exhibit scaling, which is similar to a typical "DSS scaling" . Figs. 15.251 
and 15.261 show the scaling of the scalar curvature for the families ()5.17|) and ()5.19|) 
aXr] = 0.16. The logarithm of the scalar curvature as a function of ln{p*—p) shows 
oscillations, but not enough of them in order to judge, whether these wiggles are 
superimposed on a straight line, or whether the wiggles are almost periodic. Using 
a higher numerical precision probably would yield a clearer picture. 

Unfortunately, we were not able to produce reliable results concerning the scaling 
of the black hole mass at 1] = 0.16. In all the evolutions, we have looked at, we 
find that there are two peaks in the function 2m/r, which come close to the 
threshold (0.995) towards the end of the evolution. The inner peak is afflicted 
with numerical errors, nevertheless it slows down the evolution, such that in some 
cases the outer peak cannot reach the threshold anymore. The result in some 
cases is a "broken" mass scaling. Clearly further work has to be invested here, 
before any conclusions can be drawn. 

For rj = 0.1726, where we assume the DSS solution to be the critical solution, 
although not approached as closely as at higher couplings, the scaling should be 
more conclusive. Fig. 15.271 shows the quantity ln(max7?.(M, 0)) for sub-critical 
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data of the family ()5.17j) . A straight hne with slope —27 and 7 = 0.1045 was 
fitted to these data with naked eye. In order to investigate the fine structure 
more closely, this straight line was subtracted from the data. The result is shown 
in Fig. 15.281 The difference ln(max7?.(u, 0)) — / is almost periodic, the period 

u 

being roughly A/27, where the echoing period A was taken from the directly 
constructed DSS solution. The reason for this periodicity not being exact might 
lie in the fact, that the DSS solution is not approached close enough. Again 
higher numerical precision could clarify things. 

The scaling of the black hole mass at 77 = 0.1726 is shown in Fig. 15.291 The 
"worms" displayed approximately align along a straight line with slope 7 = 
0.0965. The wiggles again are only close to periodic, as shown in Fig. 15.301 
We don't know, whether the discontinuities in the mass scaling in Fig. 15.291 stem 
from a systematic error in the measurement of the black hole mass. 
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Figure 5.21: For rj = 0.16 the intermediate asymptotics of near critical data 
(family ()5.17|1 with fixed width (7 = 1; number of grid points initially = 8049) 
are shown (dots; not every grid point is plotted). The four plots are snapshots 
at various times u, where u increases from left to right and from top to bottom. 
On the first and third plot the CSS solution is superimposed (lines), where it has 
been shifted (in In r) horizontally such that the first monotonic region agreed best 
with the evolved data. (The best fit was determined automatically by a fitting 
script by J. Thornburg). Given the horizontal shift in Inr, the corresponding 
culmination time u* is determined, as well as the location of the past SSH, which 
is denoted by a vertical fine in these plots. One clearly sees that the evolved 
data approach the CSS solution, turn away and then approach its negative. The 
culmination times associated which each CSS episode increase. 
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In(^) 

Figure 5.22: The deviation of the critical solution from the CSS solution at 
the last but one episode. Plotted is 50(r, z) = 0(r, z) — (pcss^z) at several 
time steps between r = 0.1576 and r = 1.991 (dots). As in Fig. 15.181 the 
functions a{T){yunst + b{T)y gauge) are overlaid (solid lines) with fitted values of 
the parameters a and 6, such that the maxima agree. Again the gauge mode is 
taken into account in order to correct for the uncertainty in u* . Although the 
agreement of the shapes is not very good, on can infer from Figs. 15.231 and 15.241 
that (at least for a short time) the maximum of 50 grows exponentially with a 
rate, which is close to the eigenvalue of the unstable mode A = 5.202. 
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Figure 5.23: The parameter a of Fig. 15.221 in dependence of r (dots). The sohd 
hne has a slope of 5, which is close to the eigenvalue A = 5.202 computed from 
the perturbation analysis. 



HKr)) 




Figure 5.24: The parameter b of Fig. 15.221 in dependence of r (dots). The solid 
line has a slope of 1 — 5 = 4, according to ln6(r) = (1 — A)r + const. 
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Figure 5.25: Scaling of max7^(M,0) for sub-critical data at ?7 = 0.16. The family 

u 

of initial data was ()5.17|1 with fixed width a = 1, the number of grid points was 
N = 8049. 
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Figure 5.26: The family IH". 1 91 qualitatively yields the same scaling as in Fig. 15.251 
The resolutions shown are = 500, 1000, 2000, 4000. 
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In (max 7^ (m, 0)) 
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Figure 5.27: Scaling of max7?.(M,0) for sub-critical data a.t rj = 0.1726. The 

u 

family of initial data was ()5.17|1 with fixed width a = 1, the number of grid 
points was N = 2000. 
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Figure 5.28: The straight line /(ln(p* — p)) = — 27ln(p* — p) + k was fitted to 
lnmax„7^(M, 0) in Fig. 15.271 with naked eye. The fit gave 7 ^ 0.1045. In order to 
look at the fine structure this function was subtracted from the scalar curvature. 
The result is shown in this figure ("+"). In order to check for periodicity, the 
same data were re-plotted ("x"), with a shift in ln(p* — p) of A/27, can be 
seen, the periodicity is not exact, but close. 
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Figure 5.29: Scaling of the black hole mass for super-critical data at 77 = 0.1726 
(family ()5.17p with fixed width a = 1, number of grid points was N = 2000). 
Almost all of the runs in this plot stopped because of du < 10~^^m. 




Figure 5.30: The straight line /(x) = + k with 7 = 0.0965 was subtracted 
from the data in Fig. 15.291 The result is an almost periodic function of \n{p — p*) 
with a period roughly equal to A/27. 
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Chapter 6 

Discussion and Outlook 



In this thesis we have reported on our work on the self-gravitating SU(2) a model 
in spherical symmetry. We have described our results concerning static solutions 
in the presence of a positive cosmological constant A, on self-similar solutions 
and on type II critical phenomena. 

We have shown numerically that the model (with a positive cosmological con- 
stant) admits a discrete one-parameter family of static spherically symmetric 
regular solutions. These solitonic solutions are characterized by an integer ex- 
citation number n. A given excitation will only exist up to a critical value of 
the coupling constant r]\ the higher n, the lower the corresponding critical value. 
Our calculations indicate that the infinite tower of solitons present on a de Sit- 
ter background persists at least up to a value of 77 = 1/2. Thus there exists a 
?7 > 1/2 beyond which the number of excitations is finite and decreases with the 
strength of the coupling. Qualitatively the a model under consideration shows 
striking similarities to the EYM system as studied in detail by Volkov et.al. (75j . 
The main difference being that the static solutions to the EYM-system depend 
on the value of the cosmological constant while in our case A scales out from the 
equations and t] plays the role of a "bifurcation" parameter. Another difference 
concerns the globally regular static solutions with compact spatial slices. For 
the EYM system these appear for definite values of A(?t,) while for the a model 
the corresponding solutions exist only in the (singular) limit as A goes to zero 
and definite values of rj. Thus in our case there are closed static universes with 
vanishing cosmological constant, the lowest excitation being the static Einstein 
cosmos. This is possible because in this case the stress-energy tensor of the a 
field is of the form of a perfect fluid with the equation of state p = — /i/3. An- 
other interesting aspect is the geometry of a given excitation as a function of the 
coupling strength: the static region is always surrounded by a Killing horizon 
separating the static from a dynamical region, which for small couplings becomes 
asymptotically de Sitter. As the coupling is increased the two-spheres of sym- 
metry beyond the horizon are first past and then become future trapped and 
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a cosmological singularity develops. Finally, for even stronger couplings, again 
the region beyond the horizon collapses, but within the static region the in- and 
outgoing directions (as defined by the sign of the expansion for null geodesies) 
interchange. 

An important question to be answered is whether these solitons are stable under 
small radially symmetric time dependent perturbations. Here we have presented 
preliminary results, stating that for A > all excitations are unstable with their 
number of unstable modes increasing with n. This was to be expected at least 
for small coupling. The lowest excitation thus has a single unstable mode and it 
is known, from other models, that such a solution can play the role of a critical 
solution in a full dynamical treatment of spherically symmetric collapse. 

As a further step, it would be interesting to study existence and stability of 
static solutions to this model, that have no regular center of spherical symmetry, 
but rather a static region, which is bounded by two horizons, in analogy to the 
Schwarzschild-de-Sitter (Kottler) spacetime. This is currently investigated by N. 
Miillner 

In a further part of this thesis we have numerically reproduced the results of 
Bizon and Wasserman ^1] concerning continuously self-similar solutions to the 
self-gravitating a-model in spherical symmetry. We also supplemented this work 
with a stability analysis. As in ^T] we find that a countably infinite spectrum of 
CSS solutions exists up to a maximal coupling rj = 0.5. For vanishing coupling 
T) = 0, the ground state of this family was already given in closed form by Turok 
and Spergel [73]. The homothetic Killing vector, generating continuous self- 
similarity, is timelike inside the past self-similarity horizon (the backwards light 
cone of the culmination point, where in the coupled case a space time singularity 
occurs). For small couplings all members of the family are regular up to the future 
self-similarity horizon. For larger couplings - the critical coupling depending 
on the excitation number - these solutions contain marginally trapped surfaces. 
A linear stability analysis, carried out with two different (numerical) methods, 
revealed that the number of unstable modes corresponds to the excitation number 
of the solution. In particular the ground state is stable and the first excitation 
has one unstable mode. 

The stability properties of the CSS ground state and the first CSS excitation make 
these solutions relevant for the dynamics of the system. We have shown, that for 
very small couplings the CSS ground state is the global "end state" for a set of 
initial data. The singularity at the culmination point in general is not shielded 
by a horizon, such that the CSS ground state gives rise to the formation of naked 
singularities. Although the formation of naked singularities is a general end state 
in this model for small couplings, this behavior cannot be viewed as a violation 
of the cosmic censorship hypothesis. The blow up (of the energy density) also 
occurs in fiat space, the formation of naked singularities therefore is not due to 
gravity, but to the matter model itself. On the contrary gravity regularizes these 
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singularities in the sense, that for bigger couphngs the only possible end states 
for "strong" initial data are black holes. 

The first excitation, which lies at the threshold of singularity formation in flat 
space PU] also is the critical solution between black hole formation and dispersion 
for small couplings (respectively between the formation of naked singularities and 
dispersion for very small couplings). 

For large couplings on the other hand, the solution at the threshold of black 
hole formation is discretely self-similar. Both our results of time evolution and 
critical searches on one hand and the "direct construction" using the (discrete) 
symmetry and pseudo-spectral methods on the other hand show, that the period 
A of the DSS solution rises sharply below 77 ~ 0.3. We were able to construct 
the DSS solution down to a coupling of rj = 0.1726^. At this lowest coupling 
we compared the DSS solution to the first CSS excitation, and found that there 
exists a phase of the DSS solutions, where the shapes of both functions agree 
rather well. At 77 = 0.18 the agreement is not so good. This suggests that the 
DSS solution bifurcates from the CSS solution at some coupling rjc- Due to the 
fact, that the CSS solution does not change its stability around rjc we suggest, 
that the bifurcation is a global bifurcation (in contrast to a local bifurcation), and 
furthermore that we deal with a homoclinic loop bifurcation. The theory then 
predicts, that the period A scales as the logarithm of f]c — V- fit determined 
rjc ~ 0.17. The results of the stability analysis of the DSS solution are in good 
agreement with the scaling of the black hole mass and the Ricci scalar: the DSS 
solution has one unstable mode with an eigenvalue of A ~ 9.0, which is almost 
independent of rj. At r] = 0.1726 the numerical results of the stability analysis 
are not conclusive. 

In view of the above described bifurcation scenario, we can expect the critical 
phenomena in the transition region, where the critical solution changes from CSS 
to DSS, to be rather complicated. Our results support the following view: at 
couplings ?7 > 0.17, where the DSS solution exists, it is the critical solution. Due 
to the "closeness" of the DSS and CSS solution, the CSS solution is approached 
(and left) in several "episodes" . Presumably at couplings close to rjc the stable 
modes of the DSS solution damp out more slowly than at larger couplings. With 
the given numerical precision near critical data therefore cannot approach the 
DSS solution as close as at higher couplings. This results in a periodicity of 
quantities like 2m/r, which is not exact, as well as in a scaling of the Ricci 
scalar, which has not an exactly periodic fine structure. Increasing the numerical 
precision, we could reduce the admixture of the unstable mode in the initial data, 
thereby prolonging the "life time" of the critical solution. We speculate, that then 
it would be possible to see an approach to the DSS solution which is as close as 
at larger couplings. 

^Proceeding further down would need more Fourier coefficients and therefore would increase 
the computational costs considerably. 
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At couplings below the bifurcation value of rjc ~ 0.17, we see near critical evo- 
lutions approach a configuration which still has some resemblance to discrete 
self-similarity and which shows CSS episodes. Speaking in terms of dynamical 
systems it is possible, that the critical solution now has not an exact symmetry, 
but is an invariant manifold of orbits, that "almost" close. In phase space this 
invariant manifold lies in the vicinity of the CSS solution. The smaller the cou- 
pling, the more pronounced are the CSS episodes, until at a coupling of ~ 0.14 
the CSS solution is approached only once and is the critical solution. 

To our knowledge this kind of transition from CSS to DSS in the critical solution 
and the phenomenon of a homoclinic loop bifurcation, has not been observed 
up to now in the context of type II critical collapse. Liebling and Choptuik for 
example [22] report on a transition from CSS to DSS in the Brans-Dicke modeP. 
But this transition is connected to a change in the stability of the CSS solution. 

Further work has to concretize some of the above results mainly with numerical 
methods. The code would need some improvement for the treatment of the origin, 
and it has to be tested with respect to convergence for (very) small couplings. 
Concerning the possible end states for strong initial data at very small couplings 
it would be interesting to determine the transition (in 77) from naked singularities 
to black holes more precisely. For small couplings, where the critical solution is 
CSS, it would be necessary to trace down why the critical exponent deviates by a 
few precent from the theoretically predicted value. Convergence tests and using 
a higher numerical precision could give first hints. In order to investigate the 
phenomena of episodic CSS in more detail, the first step could consist again in 
using a higher numerical precision. 

The "DSS code" could be used to study discretely self-simlar or time periodic 
solutions in other models. In particular it would be interesting to investigate, 
whether gravity is responsible for the existence of discretely self-similar solutions. 
As suggested by P. Bizon 0, a first step into this direction would be to study a 
certain artificial matter model in flat space with a self-interaction which mimics 
the interaction with gravity. 



^It might be interesting to note, that this Brans-Dicke model in spherical symmetry can be 
viewed as a ct model with a two-dimensional target manifold of constant negative curvature. 
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Appendix A 



The Shooting and Matching 
Method 



As the shooting and matching method is used at several places in this work, we 
give a short description of it here. A good description of this method can be 
found e.g. in [Hni- Consider the following ODE boundary value problem, given 
by a coupled system of ODEs, 

y\t)=F\y,t), t = l,...N. (A.l) 

and the following A^ Dirichlet boundary conditions at the ends of the interval 

[a, 6], 

gi^{y{a)) = ^, j = l,...,M, 

g,\y(h)) = Q, k = l,...,N-M. (A.2) 

So M variables at the left boundary depend on the remaining N — M variables, 
which are free parameters at, k = 1, . . . , N — M, and on the right boundary there 
are N — M variables, that depend on M free parameters bj,j = 1, . . . , M. 

After choosing the values of A^ — M variables at the left boundary and M variables 
at the right boundary, which we subsume with c = {ak,bj), Eqs. HA.l^ can be 
integrated from both ends of the interval to some matching point tmatch ^ (o, b). 
Of course the values of the solutions yitmatch) at the matching point resulting 
from the integration from left and from right will in general not agree, but there 
will be a "miss distance" 

f{c) = yieft [tmatch] a.k) " Vright {Uaatch] bj) i = 1, . . . N , k = 1, . . . N - M , 

j = l,...,M. (A.3) 

If F is smooth, then / will depend smoothly on the parameters c. The aim now 
is to find those values of the parameters for which / evaluates to zero^. This can 

^It might be convenient to replace the solutions on the right hand side of IjA.Sp by some 
function thereof 
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be achieved via a Newton iteration. For values c close to the initial guess cq, the 
miss distance can be expanded as 



If the initial guess Cq is close to a zero of / a first approximation to this zero is 
given by c with 



neglecting the higher order terms. If furthermore the Jacobian J*j(co) = 9/*/ dc^ 
is invertible we can solve for c, 



If / depends linearly on the parameters, then one step is enough, for nonlinear 
relations several steps have to be applied in order to approach the zero. 

Numerically, the Jacobian is computed by first integrating the ODEs (jA.ljl with 
the initial guess cq, and then with the perturbed values cq + (^Cfc, k = 1, . . . ,N 
and (5cfc)* = eS^k- The Jacobian can then be obtained by e.g. forward differencing 




(A.4) 




(A.5) 




(A.6) 




(A.7) 



e 
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Appendix B 

Discrete Fourier Transform 



In suitably chosen coordinates discrete self-similarity manifests itself in a peri- 
odic dependence on one of the coordinates. This suggests to work with Fourier 
expansions. The problem of constructing a DSS solution then reduces to an ODE 
boundary value problem for the Fourier coefficients. In Sec. lB.ll we review the ba- 
sic properties of Fourier series and truncated Fourier series. In practice we don't 
work with truncated Fourier series, but with the discrete Fourier transform, which 
establishes a relation between discrete Fourier coefficients and the values of the 
function at "grid points" in real space. The discrete Fourier transform can be 
viewed as the discrete approximation to the (continuous) Fourier transform, (see 
Sec. IB.2jl . In Sec. IB. 31 we define differentiation within this framework. Sec. IB. 41 
explains how algebraic manipulations are executed in "real" space, involving a 
pair of forward and backward transformations. There has to be taken special 
care in order to reduce aliasing errors, that result from such a process. 

This appendix follows closely ^H], which gives a good description of discrete 
Fourier transform, pseudo-spectral methods and aliasing errors. We only give 
the basic definitions and cite the main results. For further details and proofs we 
refer to Our discussion concerns functions defined on the interval [0, 27r]. For 
functions, that are defined on the interval [0, A] any occurrence of the independent 
variable x has to be replaced by 27rx/A. 

B.l Truncated Fourier Series 

Consider the Hilbert space of (Lebesgue-) square integrable functions L^(0,27r) 
with scalar product 




(B.l) 
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and norm 

27r 



\\u\\ = J \u{x)\' dx. (B.2) 



The functions (f)k{x) = e^^^, A; G Z form an orthogonal system with respect to 
this scalar product IB. II 

27r 

bk{x)(j)i{x) dx = 2TTSki. (B.3) 



For u G L'^{0, 27r) the Fourier coefficients of u are given by 

u^ = — f u{x)e-'^'' dx keZ. (B.4) 

277 J 



If u is real then u^k = Uk- 

The (formal) Fourier series of u is given by 

oo 

Su{x) = ^ Uk(t)k{x). (B.5) 

fe=— oo 

The truncated Fourier series of order is the trigonometric polynomial of degree 
N/2 

N/2-l 

Pnu{x) = ^ke'^'' (B-6) 

k=-N/2 

Defining the space of trigonometric polynomials of degree N/2 as 

Sn = span{e'''''\ - N/2 < k < N/2 - 1} (B.7) 

P]\fU is the orthogonal projection of u upon the space Sj\f with respect to the 
scalar product IB.H 

{Pnu,v) = {u,v) for all v G Sn. (B.8) 

Equivalently P^u is the closest approximation of u within with respect to the 
norm ESI 



For u G L^(0,27r) the Fourier series of u, Su{x) converges to u in the norm 
lR2l that is 

\u{x) - Pn{x)\^ dx asiV^oo. (B.9) 



^The convention to discuss truncated Fourier series in terms of the trigonometric polynomial 
of degree N/2 rather than N is not common in the literature but is special to ^B] 
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The Parseval identity states that 



\u\ 



|2 



2n J2 l^fcl'' (B-10) 



in particular the series on the right hand side converges. If u satisfies additional 
criteria, the convergence IB. 91 can be improved. E.g. if u is continuous, periodic 
(u(0"'") = u(27r~)), and of bounded variation on [0, 27r], then Su is uniformly 
convergent, i.e. 

max \u{x) — Pn{x)\ ^ a.s N oo. (B-H) 

a:e[0,27r] 

Concerning the rate of convergence the Parseval ident itv IB . 1 01 gives the following 

||u-P^|| = (2n \uk\^y^^. (B.12) 

k<-N/2 
k>N/2 

On the other hand for u sufficiently smooth and periodic we have 



max \u(x) — Pn(x)\ < > \uk\ ■ (B.13) 

X e [0,27r] 

k>N/2 

So the rate of convergence of the Fourier series is connected to how fast the 
Fourier coefficients of u decay. 

We are interested in the following result: if u is smooth (C°°) and periodic with 
all its derivatives on [0, 27r] then the Fourier coefficients Uk decay faster than any 
negative power of k. Of course this only applies for k bigger than some ko, the 
minimal frequency which is needed to represent the essential structure of u. 

Combining this with the formulae for the error ()B.12j) . ()B.13j) one finds that for 
u satisfying the above conditions the error between u and the truncated Fourier 
series decays faster than any negative power of A^. This is called spectral accuracy, 
exponential convergence or infinite order accuracy. 



B.2 Discrete Fourier Transform 

For any integer > consider the "grid points" 

x, = ^ j = 0,...,Ar-l, (B.14) 

where for our purposes we assume A^ to be even. If u is known at these grid 
points, then the discrete Fourier transform (DFT) of u is given by 



Uk = ^Yl -N/2<k< N/2 - 1. (B.15) 



j=0 
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As the satisfy the orthogonahty relation 

N-l 



j=0 

the inverse transform is given by 

N/2~l 



N ^ otherwise, ^ ' 



J2 ^fce'""^ j = 0,...,iV-l. (B.17) 

k=-N/2 



Form a computational point of view, the discrete Fourier transform ()B.15j) in- 
volves A^^ multiplications. It is therefore an "O (A^^) -pro cess" . Fortunately there 
exists a less expensive way to compute the A^ coefficients Uk, which is called Fast 
Fourier Transform (FFT). If A^ is an integer power of 2, then the computational 
costs for the discrete Fourier transform using FFT are only of order A^ log2 A^. A 
good description of the FFT can be found e.g. in 

The polynomial 

N/2-1 

Inu{x) = Ukc'^'' (B.18) 

k=-N/2 



is the N/2 degree trigonometric interpolant at the grid points lB.14[ i.e. lNu{xj) = 
u{xj). This polynomial is also called the discrete Fourier series of u. 

The discrete Fourier coefficients Uk can be regarded as a discrete approximation 
to the continuous Fourier coefficients Uk, in that using the trapezoidal rule to 
evaluate the integral |R4| gives Uk- 

The discrete approximation to the inner product IB. II on the space Sjy ()B.7|) is 
given by 

2n 

— J2u{x,)v{x,). (B.19) 

j=0 

Due to ()B.16j) it coincides with the inner product IB. II if ti, f G S]\f , i.e. 

{u,v)j^ = {u,v) for all M, f G Sat. (B.20) 



The interpolation operator 1^ can be regarded as an orthogonal projection oper- 
ator upon the space Sn with respect to the scalar product ()B.19|1 . as 

{Inu,v)j^ = {u,v)n for all f G Sat (B-21) 

trivially. Therefore Inu is the best approximation to u within the space Sn with 
respect to the norm \ \u\\n = {u, u) at. 
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The discrete Fourier transform of u ()B.15jl can be expressed in terms of the 
continuous Fourier coefficients of u ()B.4|) : if Su{xj) = u{xj) one obtains by using 
the relation ()B.16|) 

oo 

Uk = Uk+ ^ Uk+Nm- (B.22) 

m— — oc 

This means that the k-th mode of the trigonometric interpolant does not only 
depend on the k-th mode of u but also on all the {k + iVm)-th modes (which 
cannot be distinguished from the k-th frequency at the grid points ()B.14|) ). This 
effect is called aliasing. 



We can write 
where 



Inu = Pnu + Rnu, (B.23) 



oo / oo 



fc=— oo 



m— — oc 
m=0 



i?7vu is orthogonal to u — Pnu with respect to the scalar product ()B.1|) and 
therefore we have 

I |u - InuW^ = \\u- Pnu\\^ + WRnuW^. (B.25) 

So the error of the discrete Fourier series is always larger than the error of the 
truncated Fourier series, due to R]s[U, which is called aliasing error. Nevertheless 
it can be shown, that asymptotically the truncation errors and the interpolation 
errors decay at the same rate. 

The sequence of interpolation polynomials shows similar convergence properties 
as the sequence of truncated Fourier series. E.g. for u continuous, periodic and of 
bounded variation on [0, 27r], I^u converges uniformly to u on [0, 27r]. Concerning 
the fall off of the discrete Fourier coefficients we have e.g. for u being and 
periodic with all its derivatives: for any fixed k ^ Q and any positive N such 
that N/2 > \k\, let Uk = u\!^^ be the k-th Fourier coefficient of I^u- Then 
Eg. IB. 2 2l shows, that \u\!^'^\ decays faster than algebraically in k~^, uniformly in 
N. Using analogous arguments as in Sec. IB. If we therefore get, that the error 
between a periodic C°° function and its discrete Fourier series decays faster than 
any negative power of A^. 

In this work we use a slightly modified interpolating polynomial. First we define 
«o = j^Y^i^j) (B-26) 

j=0 

2 ^"^ 2 / ■ 

= ]vS''^''^'^'°'^^^ / = l,...iV/2-l (B.27) 

j=0 
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2 ^""^ 2 / ■ 

-5^^x,)sin(^) / = l,...iV/2-l (B.28) 

j=0 



aN/2 = ^ ^ ^(xj) cos(7rA;), (B.29) 



so 



Uk = ^{ak-ibk), -N/2 + l<k<N/2-l,k=^0 
Uq = ao, U-N/2 = a-iv/2- (B.30) 
As a_fc = Ofc and = — 6fc we can re-write the expression for u at the grid points 

N/2-1 N/2-1 ^ 

u{xj)=ao + akCos{kxj) + 6^ sin(A;Xj) + aAr/2 cos(— x-,). (B.31) 

fc=l k=l 

We now define the interpolating polynomial I^u to be 

N/2-1 N/2-1 ^ 

iNu{x) = ao+ akCos{kx)+ 6^ sin(fcx) + aAr/2 cos(— a;) (B.32) 

k=l k=l 

By definition Ij^u{x) agrees with I]^u{x) at the grid points, but differs from the 
latter in between (the difference arising solely in the N/2 frequency term). 



B.3 Differentiation 

The derivative of Inu{x) with respect to x is given by 

~ , TV TV 

{Inu)'{x)= y. ctA:(— fc) sin(/cx) + y, bkk cos{kx) — a]\f/2— sm{—x) (B.33) 

k=l k=l 

As sm{Yx) vanishes at the grid points Xj, we define the connection between the 
coefficients of the collocation derivative of u and the coefficients of I^u to be 

(ap)o = 

{ap)k = khk /c = l,...,A^/2- 1 

{hp)k = -kau k = l,...,N/2-l 

{ap)N/2 = 0, (B.34) 

where {ap)k and {bp)k denote the coefficients of the collocation derivative. Note 
that interpolation and differentiation do not commute, unless u G Sn- One can 
show that collocation differentiation is spectrally accurate. 
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B.4 Pseudo-spectral methods and Aliasing 



The way we construct the DSS solutions (see Sec. 14. 4|) our basic variables are the 
discrete Fourier coefficients (with respect to r) of the field its derivative 0' and 
the metric functions (3 and — . In order to set up the ODEs (in the spatial variable 
z) we have to compute the discrete Fourier coefficients of e.g. products of these 
functions or e.g. the sine of 0. This can be done by applying the inverse Fourier 
transform ()B.17|) to the coefficients, carrying out the algebraic manipulations and 
taking the sine in "r-space" and then transforming the result back to Fourier 
space. The overall computational scheme therefore includes operations carried 
out in Fourier space as well as operations carried out in "r-space" . Such a method 
is called pseudo spectral method. 

This transforming back and forth has to be carried out with care if one wants 
to keep the aliasing errors as small as possible. To see this consider the smooth 
and periodic functions u{x) and v{x), with their Taylor series expansions u{x) = 

oo oo 

ui^e^^^ and v{x) = ^ VkC'^^^. We denote their product by w{x) 

k=—oo k=—oo 

w{x) = u{x)v{x). (B.35) 
The Fourier coefficients of w, then are given by 

Wk = UmVn — OO < m,l < OO. (B.36) 

m+l=k 

Given the discrete Fourier coefficients Uk and v^, and the corresponding inverse 

N/2-l 

transforms Uj = ^ tt^e*'^^^ and the analogous expression for Vj, we define Wj 

k=-N/2 

to be their product in real space, 

Wj = UjVj. (B.37) 

N-l 

The discrete Fourier coefficients of Wj then are given by = ;^ Yl Wje~^^^^ and 
therefore 

Wk= UmVi+ - N/2<m,l,k< N/2-1. (B.38) 

m+l=k m+l=k±N 

We assume for a moment, that the coefficients Uk,Vk for k < —N/2,k > N/2 
are neghgible. Then we have Uk — Uk, Vk — Vk, i.e. the aliasing errors due to 
interpolation are negligible, and the sum in ()B.36j) equals the first sum in ()B.38|1 
for —N/2 < k < N/2. Then clearly the second sum in ()B.38j) introduces an error 
into the product, which again is called aliasing error. 
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In the following, we describe a method, which reduces this error. The basic idea 
is to increase the number of Fourier coefficients to M before transforming to 
real space, carry out the manipulations in real space with this higher number of 
grid points, transform back to Fourier space and then throw away the additional 
modes. If M is chosen appropriately the second sum in ()B.38j) does not contribute 
to the relevant frequencies of w. 

Let M > N. We define the following M discrete Fourier coefficients 



and the analogous expressions for V^. Given these Fourier coefficients, we get the 

M/2-1 

inverse transforms Uj = U^e^^^^ where j now runs from to M — 1 and 

k=-M/2 

the grid points Xj are given by xj = 27r j'/M and the analogous expression for 
Vj. Let again Wj = UjVj and Wk = jj Yl Wje~^''^^ , then the coefficients Wk are 
given according to ()B.38jl 



Wk= ^™^'+ Yl -M/2<m,l,k<M/2-l. (B.40) 



The strategy now is to consider only the coefficients Wk for —N/2 < k < N/2 — 1 
and throw away the higher frequencies. As Um and Vi are nonzero only for 
-N/2 <m,l< N/2 - 1, the first sum in (lR4n|l equals the first sum in (|R38|) 
for —N/2 < k < N/2 — 1. The aim is now to choose M such, that the second 
sum in flB.40j) does not contribute to the frequencies of interest. We have —N< 
m + l < N-2 and M -N/2 < k + M < M + N/2-1 and -M-N/2 < k-M < 
-M + N/2 - 1. So if M > 3A^/2 - 1 the second sum in does not contribute 

to Wk for -N/2 <k< N/2-1. 

For products of higher order M has to be chosen larger. 




(B.39) 



m,+l=k 



m+l=k±M 
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Appendix C 



The "Diamond-Integral- 
Characteristic-Evolution" Code — 
DICE 



The DICE code evolves the self-gravitating SU(2) a model in spherical symmetry 
(optionally with cosmological constant)^. Given initial data ipoir) = (j){uQ,r) at 
the initial null slice, Eqs. ()2.48|1 . ()2.55|1 and ()2.56|1 are solved numerically, while 
the "additional" Einstein equations, Eqs. ()2.57|1 and ()2.58|1 merely serve for tests 
of consistency and accuracy (see Sec. l(I5j) . The diamond integral scheme by 
Gomez and Winicour j^Sl IHH is used to integrate the wave equation ()2.48|) 
(see Sec. IC.1|) . Grid points freely fall along ingoing null geodesies. The hyper- 
surface equations ()2.55|) and ()2.56|) as well as the geodesic equation ()2.2U|) in the 
latest version are integrated using a second order iterated Runge-Kutta scheme 
(see Sec. I(I3|1 . In the vicinity of the center of spherical symmetry the integra- 
tion schemes (except for the geodesic equation) are replaced by Taylor series 
expansions (this follows [21]; see Sec. IC.2|) . 

The "kernel" of the DICE code was originally developed by Sascha Husa. Fur- 
ther improvement, development and testing of the code - including improvement 
of the integration scheme for the ODEs, measurement of the black hole mass, 
convergence tests, critical search modus etc. - was done by Michael Piirrer and 
Jonathan Thornburg. M. Piirrer also added the feature to evolve the massless 
Klein Gordon field with a compactified radial coordinate (see [Hlj). I helped to 
develop the physical and analytical foundations of the code, but did not actually 
participate in writing the code. A description of the DICE code can be found in 
US] and in 

^It also optionally evolves the self-gravitating massless Klein-Gordon field in spherical sym- 
metry 



151 



C.l The NSWE Algorithm 



The "central" evolution algorithm uses the NSWE scheme by Gomez and Wini- 
cour [SniEIlES- It is based on the fact, that the wave operator Hg in spherical 
symmetry can be expressed in terms of the wave operator in 2-dimensions Uh, 
which is defined with respect to the two dimensional metric 



dsl = -e^f^du {^du + 2rfr^ . (C.l) 



Setting if) := r0, we can write 

□,0 = -a.^-^(^9.-j^. (C.2) 

Now we use the fact, that any two dimensional metric is conformally flat, i.e. by 
setting du = ydu we have 

dsl = '^^^h with dsf^ = —du {du + 2dr) . (C.3) 

r 

The wave operator transforms under this conformal transformation as 

□^^ = e-2^^a^7/.. (C.4) 

Consider now the parallelogram S spanned by the four null lines u = Uo,u = 
Ui,v = Vo,v = vl (see Fig. IC.lp . If we integrate Eq. IC. 41 over S we get 

J d'^xV^Uh'ip = j d^x\/-hUj;tp. (C.5) 

S E 

In double null coordinates (tt, ?) = u + 2r), reads 

□^V = -^dud,i), (C.6) 

so ()(y.5|l gives 

d'^xV—hnhip = ~2 J J dudvdudvip = 

S MO So 

= 2{-^PN + i^w + i^E-iJs)- (C.7) 



is chosen such that it is bounded by two null slices, separated by one numerical time 
step, and the two ingoing null geodesies along which two neighbouring grid points move. 



152 



From Eqs. (j2^ and ()(I2|1 we have 



r \ r I T 



and therefore 



^Af = ipw + '4'E-i's-^ I d xV-h{-e ^{dr—]i' + 



^w + ^E-i^s-^ 1 1 du dr {- (dr-) ^ + ^^^^^^^^^ ] . (C.9) 



The integral on the right hand side can be approximated to second order by 

du drf{u,r) ^\Ue + fw) Ar, (C.IO) 
where Am = u^^^ — and Ar = \{rE — rs + rjsi — rw)- 

Assuming now, that the fields /3, and are known at the points S", W, 
and that furthermore the r coordinate of is known, then the field ip at N can 
be computed via lC.9l and lC.10l 



C.2 Treatment of the Origin 

The values of the field ip and the metric functions j3 and — close to the origin 
are computed using a Taylor series expansion. Using the fact that ip{u, 0) = 
0) = 0, duil>{u, 0) = duip'{u, 0) = 0, etc. and the field equation ()2.48j] . V'(^o + 
Am, r) is given by 

^K + AM,r) = l^"(«o,0)r2 + V(no,0)r3 + V(Mo,0)r2An + 0(A^). (C.ll) 
zoo 

Using the hypersurface Eqs. ()2.55|1 and ()2.5(ip together with the Taylor expansion 
(jCIllj) the expansions of (3 and ^ at the time uq + An read 

(3{uo + Am, r) = r/ ( ^V^"(mo, 0) + ^V^"(mo, 0)V^'"(mo, ^Y^du + 

+ ^^"(no, 0)^'"(no, 0)r3 ) + O(A^) (C.12) 

r V 8 12 

+ ^V^"(Mo,OX'(Mo,0)r3) +0(A^). (C.13) 

At the time step the coefficients iIj"{u'',0) and ip'"{u^^Q) are extracted from 
by fitting the cubic polynomial Cir^ + C2r^ to ip at the 5 innermost non origin 
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U = Ui, u = w 




Origin 

Figure C.l: A schematic diagram illustrating the NSWE-algorithm. The two 
dimensional wave operator in flat space acting oji ip = rcj) is integrated over 
the null parallelogram E. ipjsi therefore is given by if) at the points S, E, W plus 
an integral over E. 

grid points (a further version uses a fit at the time steps vJ' and u^~^ to compute 

7/'"(M^o) and ^/'"'(M^o)). 

After integrating the geodesic equation ()2.2m) the field at the first three non 
origin grid points at the time step u^^^ is computed from ()C.11|) . The metric 
functions (3 and ^ at these grid points then are given by (jC.12j) . 



C.3 Integrating the Hypersurface Equations and 
the Geodesic Equation 

As the hypersurface equations ()2.55|) and ()2.56|) at any fixed slice u = are 
ODEs, they can be integrated using a second order iterated Runge-Kutta scheme 
jiHj . For a general ODE ij = F{y,t), defining 

y';X\ = y' + ^tny',t% (c.i4) 
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yk + 1 Jg gjygj^ g^g follows 

y'"-' = y' + {F{y\t') + F{ylX\,t'^')) ■ (C.15) 

At the initial hypersurface, given the initial configuration 0(u = 0, r), the metric 
functions can be computed by integrating the hypersurface equations with the 
above scheme. At any further slice u^, the field ip has to be computed at the z-th 
grid point, using the NSWE scheme, before the metric functions can be updated 
at this grid point, using ()(I14jl and ()(115jl . 

The geodesic equation ()2.2()j) . again an ODE, also is integrated with the above 
method. One difficulty arises due to the fact, that at the time, we integrate the 
geodesic, the metric function ^ is not yet known at the next time step. This is 
solved by taking ^ and ^' at the grid point i — 1 at slice u^^^ and doing a linear 
extrapolation to the i-th grid point. 



C.4 Grid Refinements and Adaptive Time Steps 

In order to be able to "see" type II critical collapse, where the solution at the 
threshold of black hole formation is self-similar, it is essential for the code to 
resolve widely varying scales both in space and in time. There are two features 
that enable the DICE code to manage this resolution: first, the grid points freely 
fall along the ingoing null geodesies v = const. These geodesies tend to focus in 
regions of strong curvature, as they necessarily occur if the evolved solution stays 
close to a self-similar solution for some time. This helps to increase the spatial 
resolution in such regions. Second, as the grid points eventually hit the origin 
and are dropped from the grid, the number of grid points is doubled each time 
half of the grid points are lost (this follows jSEj) (the values of the grid functions 
at these additional grid points are obtained by interpolation). This method is 
most effective, if in addition the outer boundary of the grid is fine tuned to be 
(slightly outside but) close to the past SSH of the self-similar solution (again this 
was used by Garfinkle in 26 ). 

Concerning the time step, the scheme used is not restricted by a Courant-Friedrichs- 
Lewy (CFL) limit, as the numerical domain of dependence always equals the phys- 
ical domain of dependence of the grid. Nevertheless, we need an upper bound for 
the time steps in order to get enough resolution in time. Following Refs. jHOl 12^] 
this is achieved by requiring 

CAr 

Am < — (C.16) 

r 

for all grid points. This restricts the time steps such that no grid point is allowed 
to fall inwards more than C/2 grid spacings within a single time step. (Usually 
C is set to 1.5). 
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Apart from providing enough resolution in time this restriction on the time step 
also prevents the scheme to run into an apparent horizon G+ = 0, where the code 
would crash due to /5 — > oo. Instead the time evolution slows down before the 
formation of an apparent horizon due to both Ar going to zero at the apparent 
horizon (because r fails to parameterize the outgoing null rays u = const there) 
and an increase in ^ . This increase in ^ at a null slice close to an apparent horizon 
can be explained as follows: first of all according to Eq. 12.551 3 is non decreasing 
with r. So if it gets large somewhere on the slice it stays large for all larger 
r. Furthermore at large enough r the field and its derivative are small, and 
therefore the quantity 2m/r is decreasing. ^ can be written as = e^^(l— 2m/r), 
so increases monotonically and gets large outside the outermost peak of 2m/ r. 



C.5 Diagnostics, Accuracy and Convergence Tests 

In order to test the accuracy of the code we use the following quantities. First 
we compute the mass function 12.341 in two different ways, namely 

and 



2 V r 



^ + (C.18 



2 J \ r r 



For a solution to the Einstein equations, both expressions are identical, for a 
numerically computed solution nevertheless, they will differ by a small amount 
due to finite differencing errors. Defining 

5m = (C.19) 

'^total , initial 

where mtotai, initial is the mass contained in the grid at the initial slice, we have a 
measure of accuracy, which should always be <C 1. 

The other quantities that serve as a check for accuracy are the Einstein equations 
Euur ()2.57j) and -E'flfl fl2.58j) . which are not used to compute the solution. Here the 
question of normalization remains open in part. One possibility would be to de- 
fine the sum of the absolute values of each term contributing to the expressions 
Euur and Eqq, calling them E\uur\)E\0Q\ and use this as the normalization (this 
is e.g. done in 3+1 numerical relativity). The expressions Euur/ E\uur\, Eee/ E\g g\ 
then should always be -C 1. J. Thornburg chose a different normalization, which 
does not only consist of the sum of the absolute values but taking the maximum 
over the whole slice of the absolute value of each term, and then taking the sum. 
Calling this expression Emax\uur\ he uses Euur/{^ + Emax\uur\) and the analogous 
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expression for Eqq. The reason for taking the maxima was, that as we are working 
on null slices and integrate outwards from the origin, errors occurring at some 
location in the slice are transported outwards. The additional 1 in the denomina- 
tor reduces the error norm to the pure expression E^ut near the origin, and was 
introduced to make this norm better behaved close to the origin. What speaks 
in favor of this normalization is the fact, that the error norms constructed like 
this in most of the experienced situations are small, when the solution is well 
behaved, whereas they get large, when something is wrong. 

A number of convergence tests was done by J. Thornburg and M. Piirrer to test, 
whether the accuracy really corresponds to second order. Assuming that the 
fields are globally second order accurate, i.e. that the quantity ^^at, computed 
numerically with grid points differs from the true continuum solution \l/ by 
= \l/ + c/N"^, where c is independent of the resolution, and all higher order 
terms are neglected, then '^2N, computed with twice the resolution differs from 
\I' by \E'2Ar = \E' + c/ (4A^^)). Assume now that \E' = ^, as is the case e.g. for the 
above defined quantities 5m, Euur and E00. Then 

= ^ = l^N. (C.21) 

Three excellent examples of these convergence tests can be found in the appendix 
of where second order convergence is shown for the quantity Sm for near 
critical solutions, and for the quantity p*. Further convergence tests of the DICE 
code can be found in 



An additional test on the code was done by M. Piirrer who re-investigated critical 
collapse of the massless Klein-Gordon field with the DICE code, and found the 
critical solution to be DSS with the reported [IHI values of A and 7. The fact 
that the Klein-Gordon DSS solution can be resolved with this code is a good test 
for resolution, as the echoing period A ^ 3.43 in this case is much bigger than for 
the a model for large couplings, and therefore is harder to resolve numerically. 

We also mention the fact that the critical solution of the a-model for large 
couplings agrees rather well with the directly constructed DSS solutions (see 
Figs. 15.111 and I5.12|l . This provides a good test for the DICE code as well as for 
the "DSS code". 



■^If the value of ^ is not known in advance, one has to compare three different solutions, 
with N , 2N and AN grid points. Second order convergence then is given if 

^iN-^2N^Y^-^^li^2N-^N). (C.20) 
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C.6 Measurement of the black hole mass 



The formation of a black hole is signalled by the formation of an apparent horizon 
©+ = (or /3 ^ oo or ^ — > 1). As explained in the last section, the code slows 
down and stops before an apparent horizon forms. The code "detects" the black 
hole whenever ^ exceeds some threshold close to 1 anywhere in the grid. This 
threshold is usually set to be 0.995. At each time step the "momentary" black 
hole mass is defined as the mass thms of the outermost such grid point. At 
each further time step this "momentary mass" should increase, corresponding to 
additional matter falling into the black hole. If the numerical grid extends to 
large enough values of r initially and if the mass does not change substantially 
from one time step to the other, we have a good estimate for the "final" black 
hole mass. After detecting a black hole the code either stops because the number 
of time steps (after black hole detection) exceeds some hmit or because the time 
steps shrink to du/u < 10^^^ (which corresponds to the order of floating point 
roundoff errors). 

For large couplings the runs were made with an initial spatial extension of the 
grid, which lay substantially outside the backwards light cone of the critical 
solution. For smaller couplings on the other hand, the outer boundary of the 
grid had to be flne tuned to be close to the backwards light cone, in order to get 
enough resolution. In this case clearly the flnal numerical estimate of the mass 
does not correspond to the final mass of the black hole, but rather measures the 
apparent horizon close to the past SSH of the critical solution. 
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Appendix D 
Conventions 



Conventions concerning curvature quantities and the signature of the metric are 
the same as in [TH]. In particular we have 

Signature: 

(- + ++) (D-1) 

Riemann tensor: 

Tpcr ^ -pa _ -pa I pa p(7 _ pa p(T /p) 2] 

p^u i^p,f^ MPi*^ o^A* MP <^'^ \ I 

Ricci tensor: 
Scalar curvature: 

7^ = ^7^'^7^,, (d.4) 

Action: 



S = jd^x^{^{n-2A) + CM} 
Stress energy tensor: 



where dgf"^ = {dg-^Y" 
Einstein equations: 



(D.5) 



2 j d''x5g{^gCM) = j d^x^g{-T^,)5g^^, (D.6) 



G^y + Xg^u = i^T^u (D.7) 



Throghout this work the speed of light is set to unity, c = 1. 
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